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Preface 


The  purpose  of  this  study  was  to  modify  a  C02  laser  design  program  so  that  it  would  calculate 
pumping  terms  through  numerical  solution  of  the  collisional  Boltzmann  equation.  Both  the  laser 
design  code  and  the  Boltzmann  solver  routine  were  written  in  Microsoft  QuicicBasic  for  use  on  an 
IBM  PC  or  equivalent.  Transport  coefficients  were  computed  from  the  electron  number  density 
distributions  output  by  the  Boltzmann  routine. 

Four  methods  of  solving  the  Boltzmann  equation  were  tried,  with  L-U  decomposition  giving 
the  most  accurate  results  and  the  shortest  run  times.  Good  agreement  was  found  between  the 
predicted  transport  coefficients  and  the  results  reported  by  others.  Areas  where  the  program  had 
difficulty  were  documented  and  should  be  more  thoroughly  investigated.  A  few  of  the  effects  on 
laser  output  of  Boltzmann  related  parameters  were  also  studied. 

In  performing  this  project  I  had  a  great  deal  of  help  and  guidance  from  others.  I  wish  to 
thank  my  thesis  advisor,  Maj  Stone,  for  starting  this  project  and  for  helping  me  solve  some  very 
tough  problems.  I  also  wish  to  thank  Dr  Bailey  for  the  many  discussions  we  had  and  the  insights 
he  provided.  Finally,  I  wish  to  thank  my  wife  Cynthia  for  all  her  support  and  understanding  during 
this  project. 


David  Alan  Honey 
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Abstract 

The  collisional  Boltzmann  equaiton  was  solved  numerically  to  obtain  excitation  rates  for  use 
in  a  C02  laser  design  program.  The  program  was  written  in  Microsoft  QuickBasic  for  use  on  the 
IBM  Personal  Computer  or  equivalent. 

Program  validation  involved  comparisons  of  computed  transport  coefficients  with  experimen¬ 
tal  data  and  previous  theoretical  work.  Four  different  numerical  algorithms  were  evaluated  in 
terms  of  accuracy  and  efficiency.  L-U  decomposition  was  identified  as  the  preferred  approach.  The 
calculated  transport  coefficients  were  found  to  agree  with  empirical  data  within  one  to  five  percent. 

The  program  was  integrated  into  a  C02  laser  design  program.  Studies  were  then  performed 
to  evaluate  the  effects  on  predicted  laser  output  power  and  energy  density  as  parameters  affecting 
electron  kinetics  were  changed.  Plotting  routines  were  written  for  both  programs. 


A  Numerical  Solution  to  the  Boltzmann  Equation 
for  Use  in  Calculating  Pumping  Rates 
in  a  C02  Discharge  Liiser 


/.  The  C02  Laser  Model 

The  C02  laser  uses  the  vibrational-rotational  energy  transitions  of  the  C02  molecule.  This 
molecule  is  a  linear  triatomic  (1:28).  As  such,  it  is  subject  to  asymmetric  stretch,  symmetric  stretch 
and  bending  vibrational  modes  (1:29),  (2:25),  (3:8).  These  vibrational  motions  may  be  denoted  as 
CO2(100)  or  j/1,  CO2(010)  or  i/2,  and  CO2(001)  or  «/3. 

The  kinetics  of  C02  laser  has  been  modeled  as  a  four  level  system(l:28).  In  this  model,  the 
1/3  vibrational  mode  is  the  upper  laser  level,  the  ul  mode  is  the  lower  laser  level  (for  10.6  micron 
radiation),  and  the  j/2  mode  is  an  intermediate  level  through  which  the  C02  molecules  must  travel 
before  reaching  the  ground  state.  In  order  to  aid  excitation  of  the  upper  laser  level  i/3,  nitrogen  is 
added.  The  energy  difference  between  the  v=l  level  for  N2  and  the  v=3  level  for  C02  is  18cm~T 
This  energy  resonance  allows  for  rapid  energy  transfer  from  the  vibrationally  excited  states  of  N2 
(which  are  nearly  equally  spaced  (1:27))  to  i/=3  of  C02.  Because  of  its  cross  section,  N2  is  readily 
excited  in  an  electric  discharge,  and  it  does  not  radiatively  relax  (1:29). 

In  a  dc  electric  glow  discharge,  electrons  are  accelerated  by  the  applied  field.  This  directed 
velocity  is  thermalized  in  momentum  transfer  collisions  and  causes  the  electrons  to  excite  vibrational 
levels  of  C02  and  N2  through  impact.  The  rate  at  which  the  various  energy  levels  of  these  two 
species  are  populated  is  determined  by  the  velocity  distribution  of  the  electrons.  The  velocity 
distribution  of  the  electrons  is  affected  by  such  factors  as  the  strength  of  the  applied  field  divided  by 
the  neutral  gas  number  density  (E/N),  ambient  gas  temperature,  momentum  transfer  and  inelastic 
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processes  as  well  as  excited  state  populations.  The  Boltzmann  equation  relates  all  of  these  various 
factors  and  allows  computation  of  the  electron  number  density  as  a  function  of  velocity  (or  energy). 
The  next  section  develops  a  special  case  of  the  Boltzmann  equation  using  assumptions  appropriate 
to  a  laser.  Later  sections  will  apply  it  to  compute  pumping  terms  for  a  set  of  C02  laser  rate 
equations. 
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II.  The  Boltzmann  Equation 


Solving  the  Boltzmann  equation  allows  us  to  calculate  the  partitioning  of  energy  into  loss 
channels.  We  will  first  show  a  general  form  of  the  Boltzmann  equation  and  an  approximation  to  it 
by  way  of  a  two  term  spherical  harmonic  expansion.  We  will  then  summarize  the  development  of 
a  form  of  the  Boltzmann  equation  which  is  specific  to  our  interests  and  solve  it  numerically. 


2.1  Derivation  of  the  Boltzmann  Equation  and  Important  Assumptions 


We  begin  by  noting  that  the  number  of  particles  in  a  small  volume  dxdv  of  phase  space  is 


dn  =  f(x,  V,  t)dxdv 


(1) 


The  function  /(£,  v,  t)  is  the  density  of  the  particles  in  phase  space,  and  is  a  distribution  function. 
The  time  evolution  of  this  distribution  function  will  be  given  by 


dt  dt  ^  dx  dt  ^  dy  dt  ^  dz  dt  ^  dvg  dt  ^  dvy  dt  ^  dv,  dt 


using 


dv  _  P 
dt  m  dt 


dt 


2L 

dt 


+  vVrf  + 


L 

m 


(3) 

(4) 


When  electrons  and  molecules  collide,  scattering  will  occur.  This  will  change  the  electron  velocity 
distribution  function.  We  will  write  the  net  gain  of  electrons  of  a  small  volume  of  phase  space 
by  {^)cdxdv.  The  rate  of  change  of  the  distribution  function  f  as  n  result  of  scattering  through 
electron-molecule  collisions  is  given  by  j)-  Therefore 


df  , 


V  /e 


(5) 
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Since  we  are  assuming  the  case  of  a  spatially  uniform  plasma,  we  have 


V,/  =  0 


therefore 

f{x,v,t)  =  f{v,t)  ;  (6) 


and 


dt  m' 


(7) 


An  electric  field  exerts  a  force  F  on  the  electrons  in  a  plasma  such  that  the  acceleration  is  equal  to 
-eE/m.  The  Boltzmann  equation  now  becomes 


m 


(8) 


Under  the  conditions  (4:370)  of  a  small  applied  electric  field  ,  /  will  be  nearly  spherically 
symmetric  in  velocity  space.  Since  the  difference  between  the  masses  of  an  electron  and  molecule 
is  very  large,  an  elastic  collision  between  the  two  will  result  in  a  small  amount  of  electron  energy 
being  transferred  and  a  large  deflection  for  the  electron.  We  will  assume  that  the  electron  mean  free 
path  is  small  compared  to  the  physical  dimensions  of  the  laser  tube.  For  a  typical  case  of  a  laser 
tube  with  the  neutral  gas  at  standard  temperature  and  pressure,  and  assuming  an  average  cross 
section  of  5  x  10“**cm’,  the  mean  free  path  would  be  l/{NmQm)  which  is  equal  to  7  x  10~®cm. 
Most  lasers  would  greatly  exceed  this  dimension.  The  electron  velocity  distribution  will  then  be 
nearly  spherically  symmetric.  The  distribution  function  may  then  be  approximated  (5:27)  by  the 
first  two  terms  of  a  spherical  harmonic  expansion 


(9) 
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It  will  be  assumed  that  |/i|  <C  fo- 


The  collision  frequency  for  elastic  electron-molecule  collisions  is  given  by 


l/c{v)  =  NQmV 


(10) 


where  N  is  the  neutral  gas  density  and  Qm  is  the  momentum  transfer  cross  section.  Using  (9)  and 
(10),  equation  (8)  then  becomes  two  coupled  equations  (3:144) 


dfi  eSdfo  M  n  (  \r 


_ e_^ 

dt  3mt;^  dv 


(v^E-h) 


•{%). 


(11) 


(12) 


Solving  this  coupled  pair  of  differential  equations  forms  the  basis  for  our  solution  to  the  Boltzmann 
equation. 


2.2  Energy  Transfer  by  Elastic  Collisions,  the  Applied  Electric  Fields,  and  Inelastic  Collisions 

2.2.1  .Elastic  Collisions  The  rate  at  which  energy  will  flow  from  the  electrons  to  the 
molecules  as  a  result  of  elastic  collisions  is  given  by  (2:85,3:145)  the  energy  lost  by  the  electron 
(AE)  multiplied  by  the  elastic  collisions  cross  section  I{0,  |v  -  K|),  the  neutral  gas  density  N,  the 
electron  number  density  rio,  the  electron-molecule  relative  velocity  jiT  —  V\,  the  molecular  velocity 
distribution  /m(^)  the  electron  velocity  distribution  fo{v),  ^1  of  which  is  then  integrated  over 
the  electron  and  molecular  velocities  (v  and  V).  The  result  will  be  (2:88,  3:146) 


Qm{v)NnofM{V) 


fo{v)  + 


V  dv 


v^dvd^ 


(13) 


where  m  is  the  electron  mass,  M  is  the  molecular  mass,  the  relative  velocity  v  has  been  used  for 
if-  V',  V  is  assumed  to  be  much  greater  than  V. 
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The  molecuiar  velocity  distribution  function  will  be  assumed  to  be  Maxwellian  and  normalized 
to  unity 

j  MV)(^V  =  1  (14) 

The  rate  of  energy  transfer  for  elastic  collisions  then  becomes 

where  T  is  the  ambient  gas  temperature. 

£.2.2  Applied  Electric  Field  The  applied  electric  field  is  the  mechanism  by  which  the  elec¬ 
trons  in  the  plasma  are  accelerated  and  the  excitation  process  is  initiated.  An  electron  with  velocity 
V  will  have  a  kinetic  energy  given  by 

e-^rnv^  (16) 

The  rate  of  change  of  the  kinetic  energy  is 

^  =  miP  v  = —e£  ■  V  (17) 

di 

If  we  now  examine  the  total  rate  of  change  of  the  electrons  over  all  electron  velocities  due  to 
the  applied  field,  we  see  that 

W,=-jnoeS-  ^  .  A(r)]  dv  (18) 

Noting  that  /a(v)  is  an  isotropic  function  and  using  the  approximation  (2:81) 

(19) 
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results  in 


J eE  ■  fiv^dv 


(20) 


Once  steady  state  has  been  reached,  the  power  transfer  from  the  field  to  the  electrons  will 
equal  power  transfer  from  the  field  to  the  molecules  (through  elastic  and  inelastic  collisions).  We 
can  than  obtain  power  transfer  through  inelastic  collisions  only  by  subtracting  equation  (15)  from 
(20).  If  we  examine  this  result  for  the  case  of  power  transfer  through  inelastic  collisions  by  a  single 
electron  in  the  velocity  space  element  v^dv  we  obtain  the  relation 


■  h  -  ^foNQmV^  -  !^kTNQm^V^ 


(21) 


which  may  be  substituted  into  (12)  to  account  for  the  contribution  due  to  elastic  collisions.  The 
result  is 


£4 

dt 


1  d 
mv^  dv 


m 


mkT 


■j^NQ,nV*f,  +  —NQ„v^ 


dv  J  \m  J 


in€l 


(22) 


Let  us  return  for  a  moment  to  the  first  of  our  two  coupled  differential  equations 


dt  m  dv 


!^vQr„{v)fi 


A  solution  to  this  equation  (2:90,  3.T44)  is 


/;  = 


jSit 

mNvQt 


(23) 


Examination  of  the  exponential  term  for  practical  laser  plasma  parameters  shows  that  it  can  be 
neglected,  since  its  relaxation  time  would  be  very  fast  compared  to  the  time  scale  of  the  other 
processes.  For  example,  at  standard  temperature  and  pressure  we  have  (6:17)  N  =  2.69xl0'*cm“® 
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(Loschmidt’s  number).  An  electron  at  273  K  would  have  a  thermal  velocity  v  =  6.92xl08cm/sec. 
A  typical  cross  section  (6:39)  would  be  5x10" ‘*cm*.  This  would  yield  a  relaxation  time 

r  =  {NvQm)~^  -  10“*“*  seconds  (24) 


We  may  then  write 


h  = 


mNvQm 


(25) 


Using  equation  (25)  to  eliminate  /i  from  equation  (22),  and  employing  the  variable  substitution 


tl  Si 


(26) 


equation  (22)  becomes 

f2e\^d(E^u  df,  2mNQ„u^fo  2mkTNQ,nU^  dfA  _  (dfA 

dt  Vmuy  du  \3NQm  ^  du )  \  dt  J 

2.2.3  Inelastic  Collisions  The  variable  u  of  equation  (26)  is  expressed  in  units  of  volts.  I 
is  related  to  the  translational  energy  of  an  election  through 


(28) 


We  may  also  express  the  energy,  Ej,  of  the  jth  internal  state  of  a  molecule  by 


Ej  =  euj 


(29) 


The  molecular  velocities  are  assumed  to  be  very  small  in  the  system  we  are  examining, 
compared  to  the  electron  velocities.  We  will  therefore  assume  that  the  inelastic  cross  section  for 
the  excitation  of  the  molecule  from  its  ground  state  to  any  internal  state  does  not  depend  on  the 
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molecular  velocity.  The  inelastic  cross  section  for  this  j"*  internal  state,  which  we  shall  denote 
<5j(u*),  will  thus  only  depend  on  the  electron’s  initial  energy  u*  (as  determined  by  the  electron’s 
initial  velocity  u*  ). 

As  electrons  suffer  inelastic  collisions  with  and  transfer  energy  to  the  molecules,  the  distri¬ 
bution  function  at  energy  u  will  be  affected  two  ways.  First,  the  density  per  unit  volume  per  unit 
velocity  space  at  a  particular  energy  u  will  gain  as  electrons  with  an  initial  energy  of  u  +  Uj  collide 
with  a  molecule  and  transition  to  a  final  energy  of  u.  This  gain  is  given  by  (3:149) 

(^)  =  +  +  “;){«  +  «>)  (30) 

Second,  there  is  a  loss  at  energy  u  as  electrons  with  initial  energy  u  suffer  an  inelastic  collision  with 
a  molecule  and  transition  to  a  lower  final  energy  value.  This  loss  is  given  by  (3:149) 

The  rate  of  change  of  fo(u)  due  to  inelastic  processes  will  thus  be  found  by  subtracting  equation 
(30)  from  equation  (29)  to  yield 

(^)  ^  =  +  )<?>(«  + «;)(«  +  «>)- /.(w)Q;(«H  (32) 

Electrons  colliding  with  excited  molecules  may  also  gain  energy.  These  collisions  of  the  second 
kind  in  which  excitation  energy  goes  back  to  electron  are  called  superelastic,  and  the  effect  on  the 
distribution  function  is  given  by  (3:149) 

(^)  =  11  -  «;)<?->(«  -  «i)(w  -  W;)  -  /<.(")'?-;(«)«)  (33) 

where  denotes  the  cross  section  for  a  superelastic  collision. 
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We  may  now  substitute  into  equation  (27)  for  the  term  j  equations  (32)  and  (33). 


The  result  of  this  is 


3  \yV/  dxi  \Qm  du )  M  du 


Jo)  + 


2mkT  d 


Me  du 


L  ^  f 

+  ~^'^Nj[(u  -  Uj)f(u  -  Uj)Q.jiu  -  uj)  -  ufo{u)Q.j{u)]  =  )  ~ 


(34) 


This  is  the  Boltzmann  equation  for  the  isotropic  part  /»(u)  of  the  nearly  isotropic  electron  energy 
distribution  function,  for  a  gas  containing  one  type  of  molecule.  When  multiple  species  are  present, 


2mkT 
e  du 


*  J 


*  1 


1  is  (3:150) 

f] 

-  ‘a/* 

du  J 

-  uMu)Q^iu)] 

(ir)‘( 

dfo\ 
dt  ) 

(35) 

'jt/N,  and  Njt 

is  the 

number  density  of  the  excited  molecules  in  state  j  of  species  *.  The  energy  of  a  molecule  of  species 
k  in  state  j  is  given  by  uy*. 


2.3  Finite  Difference  Form  of  ike  Boltzmann  Equation 

A  numerical  solution  to  the  Boltzmann  equation  based  upon  a  finite  difference  approach  has 
been  developed  by  Rockwood  (7:2348,  8:377)  and  expanded  upon  by  Davies  (9:369).  We  will  now 
express  equation  (35)  in  finite  difference  form,  and  later  use  this  result  as  the  basis  for  our  main 
computer  program.  This  approach  consists  of  two  main  parts.  First,  equation  (35)  is  expressed  in 
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terms  of  the  normalized  electron  number  density 


n(ti)  =  u  V(«) 

/  uif(u)du  =  1 

Jo 


(36) 


Also,  a  variable  change  from  u(volts)  to  e  (electron  volts)  is  accomplished  where 


e  =  eu 


(37) 


Second,  equation  (35)  is  converted  in  to  a  set  of  K  -  coupled  ordinary  differential  equations.  This  is 
accomplished  by  dividing  the  electron  energy  axis  into  K  cells  of  width  Ac.  Having  finite  differenced 
the  energy  axis  in  to  K  cells,  we  now  project  equation  (35)  on  to  it.  The  result  is  a  set  of  K  coupled 
first  order  differential  equations,  which  have  been  approximated  by  a  finite  difference  representation. 
The  approximations  made  to  do  this  are 


dn  nt+i  -  n* 

—  T 

d(  Ae 


d^n  nt_i  -2nk  +  nk+i 

d(  ^  (Af)* 


(38) 


where  n*  is  the  electron  number  density  in  the  energy  range  {k  —  1)A«  =  f*  to  tAc^ .  The 
Boltzmann  equation  thus  becomes 


Ok-in*-!  +  bk+lTlk+l  —  (a*  +  6»)n*  +  ~  Rji.k^k) 

j  • 

+  ^  ^  "  (39) 

}  • 


where: 


2kT^\ 
Ar  j 


11 


=  (^)*2mjv5]  [i.Q;;,(.+)/A/.] 

«i..*  =  fi;.K^)  =  QJ(^?)(^) 

=  Q‘-A4Htt)  =  (^) 


+  \  i 
i 


mj,  =  the  nearest  integer  to 


eui. 


Ae 


Aj,  =  exp  {-Cj./kTJ,) 
6 


In  this  analysis: 

k  =  8.6174xio-SeV/K 
T  =  degrees  Kelvin 

Tj',=  effective  vibrational  temp  of  excited  state  j  of  species  s 
ej^  =  value  of  the  right  hand  edge  of  bin  k  in  electron  volts 
m,  Ms  are  in  kg 
Q*„,  Qj  are  in  square  meters 

Aj,  =  Nj,fN,  is  the  fraction  of  a  species  s  that  is  in  the 
excited  (internal)  state 

U;,  is  the  energy  loss  in  volts  associated  with  the  /^inelastic  process  in  species  s. 

Recalling  equations  (27)  -  (35),  the  grouping  of  terms  are  such  that  (9:98) 

at  is  the  rate  at  which  electrons  in  energy  bin  k  transition  to  energy  bin  (k+1)  due  to  elastic  collisions 
hk  is  the  rate  at  which  electrons  in  energy  bin  k  transition  to  energy  bin  (k-1)  due  to  elastic  collisions 
Ri$,k+mj.,  is  the  rate  at  which  electrons  in  energy  bin  (k  +  m^,)  transition  to  energy  bin  k 
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due  to  inelastic  collisions 


t-mj.'  which  electrons  in  energy  bin  (<r  -  my,)  transition  to  energy  bin  k 

due  to  superelastic  collisions 

VVe  will  now  apply  boundary  conditions  to  this  model.  The  bin  representing  the  lowest  energy 
will  have  index  ifc  =  1.  Since  no  electrons  can  transition  to  or  from  an  energy  bin  below  zero  energy, 
we  must  have  61  =  0  amd  Rfj,  =  0  for  t  —  mjs  <  1.  The  bin  representing  the  highest  energy  will 
have  index  k=K.  Since  no  electrons  can  transition  to  or  from  am  energy  bin  above  the  maximum 
energy  level  (e^ ),  we  must  have  a/(  =  0  and  Rj,  =  0  for  <r  +  mja  >K.  These  boundary  conditions 
were  accounted  for  in  the  computer  solution  to  equation  (39). 

2-4  Numerical  Solution  of  the  Boltzmann  Equation 

Four  numerical  were  used  to  solve  equation  (39).  These  were:  :  Gauss-Jordan,  L-U  Decom¬ 
position,  Gauss-Seidel  and  Successive  Over  Relaxation.  An  explanation  of  each  and  a  comparison 
of  the  of  these  methods  is  in  Appendix  E.  L-U  decomposition  was  found  to  give  the  best  results  for 
both  speed  and  accuracy.  The  mathematical  basis  for  these  solutions  is  given  here. 

Equation  (39)  can  be  condensed  to  matrix  form  as 


n  =  ^Cnn/  (40) 

I 

As  shown  earlier,  the  matrix  elements  of  C*/  give  the  rate  of  transfer  of  electrons  from  one  energy 
bin  to  another.  It  is  assumed  (7:2349,  8:378)  that  these  matrix  elements  are  constant.  As  can  be 
seen  from  (39),  the  matrix  C  will  be  a  banded  matrix.  The  tridiagonal  component  is  due  to  elastic 
scattering  and  the  applied  electric  field.  Elements  above  the  principle  diagonal  are  due  to  inelastic 
scattering,  those  below  are  due  to  superelastic  scattering. 
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Equation  (40)  expressed  as  a  forward  finite  difference  approximation  becomes 


nk{t  +  h)  ~  nk(t)  =  ^  hCnni(t  +  h)  (41) 

t 

where  /i  is  a  time  step.  Since  we  are  evaluating  the  right  hand  side  of  equation  (41)  at  time  (t  +  h), 
we  have  an  implicit  algorithm.  The  advantage  of  an  implicit  algorithm  over  an  explicit  one  is  that 
the  former  is  unconditionally  stable  (10:575,  11:379).  Equation  (41)  may  be  rewritten  such  that 

nk(t  +  A)  -  ^  hCnni(t  +  A)  =  nt(t) 

I 

^(/ -  hC)n,{t  +  h)  =  ni(t) 

k 

n,{t  +  h)  =  Yl(^-hC)-^nkit)  (42) 

k 

where  1  is  the  identity  matrix.  This  gives  us  an  iterative  algorithm  whereby  we  may  compute  the 
electron  number  density  distribution.  Accordingly,  we  first  form  the  C  matrix.  Next,  we  multiply 
the  C  matrix  by  the  time  step  h,  and  subtract  this  product  from  the  identity  matrix.  The  inverse 
of  this  is  then  c^dculated  and  multiplied  by  the  last  iteration  of  rik  in  order  to  determine  the  next 
estimate.  As  such  we  use  nj  to  determine  n^,  then  to  determine  nj  and  so  forth.  This  algorithm 
was  used  in  both  the  Gauss-Jordan  and  L-U  decomposition  methods. 

S.S  Transport  Coefficients 

In  order  to  provide  the  required  input  to  the  laser  design  program  C020SC,  we  must  c^dculate 
the  electron  drift  velocity  and  the  excitation  rate  of  certain  excited  states  in  specific  species.  The 
drift  velocity  is  used  to  determine  the  electron  number  density  Ne  through  (3:167) 

vseN,  =  j  (43) 
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where  j  is  the  current  density  in  the  discharge  tube.  The  product  of  the  electron  number  density 
Ne,  the  inelastic  excitation  rate  x  and  the  applicable  species  number  density  (Nco2  or 
provides  the  desired  pumping  term  for  C020SC. 


2.5.1  Drift  Velocity  Referring  to  equation  (39)  we  see  that  the  first  term  in  each  of  the 
expressions  for  at  and  bt+i  determine  the  rate  at  which  the  electrons  will  exchange  energy  with 
the  dc  electric  field.  We  will  define  new  terms  8*  and  It  such  that 


2Ne  1 

3m  ' 

St  =  ^i 

/Ey  1 

3m 

[nJ  utfN 

(44) 


The  rate  at  which  electrons  gain  energy  from  the  field  will  be 


E,  =  ^(a*  -  ht)nt^(  (45) 

k 

This  ohmic  loss  or  Joule  heating  may  be  written  as 


Eg  —  ^^(a*  —  ht)f*k^(  ^  J  .  S  —  J E 

k 


(46) 


The  units  for  this  expression  are  power  per  unit  volume  or  joules  per  cubic  meter  per  second. 
Conversion  to  units  of  electron  volts  gives  for  the  drift  velocity 


Vi 


Ylti^k  -  h)nk^( 
NeE 


(47) 


2.5.2  Rate  of  Excitation  The  excitation  rate  for  inelastic  processes  is  given  by  (7:2350) 


rji  _  f<r.,(<)v(e)n(e)d< 

'  /«(«)* 


(48) 


where  (T,j  is  the  cross  section  for  inelastic  excitation  of  the  internal  state  of  species  s.  Conversion 
to  a  formula  for  numerical  calculation  results  in 


The  units  for  this  excitation  rate  are  m^/s. 
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III.  Computer  Programs 


The  main  goal  of  this  project  was  to  write  a  computer  subprogram  which  would  enable  the 
program  C020SC  to  generate  pumping  rates  for  use  in  its  laser  kinetic  equations.  The  name  of  this 
subprogram  is  Boltz.  Both  C020SC  and  Boltz  are  written  in  Microsoft  QuickBasic  for  use  on  an 
IBM  personal  computer  or  equivalent.  QuickBasic  automatically  detects  the  presence  of  and  uses 
a  math  co-processor  (if  installed).  All  run  times  and  performance  specifications  reported  are  for  an 
Apple  IIGS  computer  with  an  Applied  Engineering  PC  Transporter  installed.  This  is  equivalent 
to  an  8  MHz  8086  computer  with  640K  of  random  access  memory  and  an  8087  math  co-processor 
installed. 

The  first  part  of  this  section  explains  the  subprogram  Boltz  and  the  second  all  of  the  subpro¬ 
grams  it  uses.  The  third  part  describes  the  cross  section  data  file  structure  and  contains  instructions 
for  user  modification  of  the  data.  The  fourth  part  lists  and  explains  modifications  to  the  program 
C020SC  that  were  necessary  in  order  to  incorporate  Boltz.  The  final  part  of  this  section  is  a  guide 
to  the  use  of  C020SC  with  Boltz. 

3.1  Boltzmann  Routine 

The  purpose  of  the  Boltz  subprogram  is  to  compute  excitation  rates  which  are  then  used  to 
calculate  the  pumping  terms  for  C020SC.  Boltz  does  this  by  first  calculating  the  electron  number 
density  distribution,  and  then  the  desired  rate  using  equation  (62).  A  flow  chart  of  this  subprogram 
is  in  Figure  1.  Our  discussion  of  Boltz  follows  this  flow  chart.  A  copy  of  Boltz  is  in  Appendix  A. 

3.1.1  Read  Parameters  Passed  from  C020SC  The  first  ten  parameters  in  the  argument  of 
subprogram  Boltz  are  all  information  supplied  to  Boltz  from  C020SC.  The  last  three  parameters 

in  the  argument  of  Boltz  are  the  output  Boltz  supplies  to  C020SC.  These  parameters  are: 

kc%  -  the  number  of  cells  or  bins  into  which  the  energy  axis  has  been  divided  as  a  result 
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of  finite  differencing 

de  -  width  of  a  bin  in  electron  volts 

pres  -  total  pressure  of  the  gas  mixture  in  atmospheres 

d2  -  percentage  of  C02  in  the  gas  mix 

dl  -  percentage  of  N2  in  the  gas  mix 

d3  -  percentage  of  He  in  the  gas  mix 

T  -  ambient  temperature  of  the  gas  mix 

VT()  -  a  three  element  vector  containing  the  vibrational  temperatures  for  the  (010), 

(100)  and  (001)  states  of  C02 

EN  -  E/N,  the  electric  field  strength  in  volts/cm  divided  by  the  gas  number  density  in  cm  , 
resulting  in  units  for  E/N  of  volts  cm“^ 

jd  -  current  density  in  the  plasma  tube  in  amps/cm^ 

Wa  -  pumping  rate  for  the  C02  (001)  level  in  r  "•^s"’ 

Wb  -  pumping  rate  for  the  C02  (100)  level  in  m“®s“' 

Wc  -  weighted  total  pumping  rate  for  ail  N2  levels  (v  =  1  to  8)  in  m~®s~^ 

3.1.2  Load  Cross  Section  Data  Before  cross  section  data  is  read  from  disk,  two  events  occur. 
First,  Boltz  checks  to  see  if  it  has  already  been  run  before.  If  Boltz  has  already  been  previously  run, 
then  the  cross  section  data  is  already  in  storage  and  need  not  be  read  again  from  disk.  This  saves 
about  seven  seconds  of  run  time.  The  variable  “event”  is  the  means  by  which  Boltz  remembers  if 
it  has  been  previously  run.  By  declaring  *ne  variable  “event”  with  the  command  STATIC,  event 
is  not  reinitialized  each  ti.„  .  is  executed.  Any  variable  declared  by  using  STATIC  will  retain 
its  value  from  the  previous  run  of  the  subprogr2un  in  which  the  variable  was  declared  and  assigned 
a  value.  Before  proceeding  any  further,  Boltz  examines  the  value  of  “event".  If  “event”  has  a  value 
not  equal  to  one,  then  the  subprogram  Boltz  has  not  been  previously  run  and  the  program  proceeds 
with  the  next  step,  which  is  to  set  “event”  equal  to  one.  If  when  evaluated  the  variable  “event” 
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does  have  a  value  equal  to  one,  then  this  is  a  signal  to  Boltz  to  skip  reading  in  the  cross  section 
data  and  instead  to  proceed  by  branching  to  the  program  marker  “label4;”.  The  second  action 
Boltz  will  perform  is  to  set  up  storage  space  for  the  cross  section  data.  These  arrays  and  variables 
are  also  declared  with  the  STATIC  command  so  that  this  data,  once  read  in  from  disk  will  not  be 
automatically  cleared  from  memory.  While  the  STATIC  command  prevents  automatic  er8tsure  of 
specified  variables  and  arrays,  it  does  not  actually  allocate  any  storage  space.  This  is  accomplished 
by  the  DIM  statement. 

The  arrays  and  variables  used  with  the  cross  section  data  are: 

eV(I,J,K)  -  energy  in  electron  volts  for  a  particular  inelastic  cross  section  K,  of  a  species  I, 
for  internal  energy  state  J 

QNj(I,J,K)  -  inelastic  cross  section  in  cm2  for  a  particular  energy  K,  of  a  species  I,  for  an 
internal  energy  state  J  (where  K  =  1  for  N2  and  K  =  2  for  C02) 
jin(I)  -  the  number  of  internal  states  for  a  species  K 

K(I,J)  -  the  number  of  pieces  of  cross  section  data  in  the  data  file  for  a  particular  species  I 
for  an  internal  state  J 

eVml(K),  eVm2(K),  eVm3(K)  -  energy  in  electron  volts  for  a  particular  momentum  transfer 
cross  section  data  specified  by  the  index  K  for  N2,  C02  and  He  respectively 

QNml(K),  QNm2(K),  QNm3(K)  -  momentum  transfer  cross  section  data  for  N2,  C02  and 
He  in  cm2,  associated  with  the  same  index  K  as  used  for  the  eVm(K)  energy  arrays 

ujs(I,J)  -  the  inelastic  threshold  energy  in  electron  volts  of  species  I  for  an  internal  state  J 
jml,  jm2,  jm3  -  the  number  of  momentum  transfer  cross  section  data  contained  in  the  data 
files  of  N2,  C02  and  He  respectively 

Other  variables  used  during  the  process  of  reading  in  cross  section  data,  but  which  are  not  static  are: 

T$  -  used  to  store  the  name  of  a  gas  as  read  in  from  the  inelastic  cross  section  data  file 
for  that  gas 
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MTS  -  used  to  store  the  name  of  a  gas  as  read  in  from  the  momentum  transfer  cross  section 
data  file  for  that  gas 

The  actual  process  of  reading  in  the  cross  section  data  file  information  for  N2  and  C02  is 
essentially  the  same.  After  associating  each  of  the  data  files  with  a  device  number  in  the  OPEN 
statements,  the  first  line  of  input  is  read  from  the  data  file.  This  first  line  tells  Boltz  how  many 
inelakstic  processes  are  in  the  data  file.  Boltz  stores  this  information  in  the  jin()  array,  and  uses  it 
in  the  first  loop  of  the  data  input  section  of  the  program  to  ensure  that  all  of  the  blocks  of  inelastic 
cross  section  data  are  read  from  the  data  file  for  that  gas.  Within  each  block  of  data  for  each 
inelastic  process,  Boltz  will  read  from  the  data  file  and  store  in  ujs()  the  inelastic  threshold  energy 
for  a  particular  internal  energy  state  and  store  in  K()  the  number  of  data  points  contained  in  the 
data  file  for  a  particular  block  of  inelastic  cross  section  data.  The  value  for  K()  is  then  used  to 
control  the  nested  loop  to  ensure  that  all  the  data  points  for  each  block  of  inelastic  cross  section 
data  associated  with  each  inelastic  process  are  read  from  the  data  file.  For  each  inelastic  cross 
section  data  point  there  is  an  energy  for  that  point,  which  is  stored  in  array  eV().  The  energy  for 
each  of  these  points  is  in  electron  volts.  The  cross  section  for  each  of  these  points  is  stored  in  th^ 
QNj()  array.  The  cross  sections,  as  read  from  the  data  file  are  in  units  of  cm*.  These  are  divided 
by  10^  to  convert  to  m*  so  as  to  be  consistent  with  the  rest  of  the  calculations. 

Momentum  transfer  cross  section  data  for  N2  and  C02  is  read  in  from  its  respective  data  file 
immediately  after  the  inelastic  data  for  that  gas.  For  He,  momentum  transfer  cross  sections  are  the 
only  data  read  in  from  disk.  The  basic  procedure  for  reading  the  momentum  transfer  cross  section 
data  for  all  three  gases  is  essentially  the  same.  The  number  of  applicable  data  points  is  read  in 
from  the  first  line  of  the  appropriate  section  of  the  data  file  and  stored  in  jml  for  N2,  jm2  for  C02, 
and  jm3  for  He.  This  value  is  then  used  to  control  the  momentum  transfer  cross  section  data  input 
loop.  The  energy  in  electron  volts  associated  with  each  data  point  is  stored  in  the  arrays  eVml(), 
eVm2(),  and  eVm3()  for  N2,  C02  and  He  respectively.  The  momentum  transfer  cross  section  data 
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are  divided  by  10^  so  as  to  convert  them  from  units  of  cm^  to  m^.  This  ends  all  transfer  of  data 
from  the  files  on  disk  to  Boltz. 

The  vibrational  temperatures  of  modes  (100),  (010)  and  (001)  are  passed  to  Boltz  from  the 
program  C020SC.  These  value  are  transferred  to  the  VB()  array  for  use  in  later  calculations. 
All  other  C02  vibrational  modes  are  set  equal  to  the  ambient  temperature.  Before  the  elements 
of  the  C-matrix  of  equation  (40)  are  calculated,  required  storage  space  is  set  aside  and  constants 
are  defined.  The  storage  space  is  declared  using  the  REDIM  statement.  Arrays  thus  defined  are 
dynamic.  Dynamic  arrays  may  be  dimensioned  using  a  variable.  For  Boltz,  this  user  defined 
variable  is  kc%,  passed  to  Boltz  from  the  main  program.  There  are  two  properties  of  the  dynamic 
arrays  used  here.  First,  it  allows  the  storage  space  to  be  tailored  to  a  particular  run  so  that 
unnecessary  storage  space  isn’t  tied  up  when  it  might  be  needed  elsewhere.  Second,  the  use  of 
dynamic  arrays  allows  for  very  large  arrays  when  compiled  with  the  /AH  option.  Under  normal 
conditions  Microsoft  QuickBasic  is  limited  in  storage  space  to  64K  per  array.  However,  the  use  the 
/AH  option  (required  for  dynamic  arrays)  allows  a  QuickBasic  program  to  have  arrays  up  to  the 
limits  of  memory  available.  The  ability  to  exceed  64K  of  random  access  memory  per  array  was 
needed  for  some  of  the  larger  arrays  reported  in  the  results  section.  For  an  uncompiled  version 
of  C020SC  run  from  the  QuickBasic  shell,  the  /AH  option  is  available  if  QuickBasic  is  started 
by  typing  QB/AH.  Versions  of  C020SC  that  were  compiled  under  a  QuickBasic  shell  that  was 
booted  with  the  /AH  option  specified,  automatically  have  this  feature  built  in.  No  special  actions 
are  required  to  run  such  a  compiled  version. 

The  arrays  defined  here  are: 

abar(I),  bbar(I)  -  these  correspond  to  values  of  a  and  b  of  equation  (56)  and  are  used  to 
calculate  the  drift  velocity  in  equation  (59)  and  the  index  I  refers  to  a  particular  bin 

a(I),  b(I)  -  these  are  the  values  of  a  and  b  from  equation  (39)  where  the  index  I  refers  to  a 
particular  energy  bin 
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Rjs(S,J.I)  -  the  rate  at  which  electrons  in  energy  bin  (k+mjs)  flow  down  to  bin  k 
due  to  inelastic  collisions  as  detailed  in  equation  (39),  and  index  S  represents  a  particular  species 
(1  =  N2,  2  =  C02),  the  index  J  corresponds  to  a  particular  internal  energy  state, 
and  the  I  index  denotes  the  energy  bin  under  consideration 

RTjs(I)  -  this  array  stores  the  sum  of  ail  the  rates  Rjs  over  all  the  species  and  all  the  internal 
energy  states  for  each  energy  bin  I,  as  dictated  by  equation  (39) 

RSjs(S,J,I)  -  the  rate  at  which  electrons  in  energy  bin  (k-mjs)  flow  up  to  energy  bin  k  as  a 
result  of  superelastic  collisions  as  shown  in  equation  (39), 
and  S,  J,  and  I  are  as  defined  above  for  Rjs() 

RTSjs(I)  -  the  sum  of  all  the  superelastic  rates  over  all  the  species  and  all  the  internal  energy 
states  for  each  energy  bin  as  indexed  by  I 

C(I,I)  -  contains  all  the  elements  of  the  C-matrix  of  equation  (40) 

NO(I)  -  the  total  electron  number  density  in  the  range  (/  -  l)Af  to  /Ae  as  shown  in 
equations  (40  -  (42) 

NOl(I)  -  temporary  storage  of  NO()  values  during  the  iteration  part  of  the  program 
PNO()  -  temporary  single  precision  storage  (for  plotting)  of  the  NO()  values,  since  Plotl  is 
single  precision,  and  a  call  to  that  subprogram  with  a  double  precision  parameter  in  the  argument 
would  result  in  a  ”  Parameter  Type  Mismatch”  error  statement 

INdx(),  vv(),  yb()  -  all  are  used  as  temporary  storage  space  in  the  L-U  decomposition  and 
backsubstitution  routines  (10:31) 

Constants  defined  in  this  section  of  the  program  are  used  to  convert  user  supplied  physical 
parameters  to  S.I.  units  and  to  supply  values  in  symbolic  notation  for  later  use  in  calculations.  The 
source  for  conversion  factors  is  reference  6.  This  part  of  the  program  also  converts  the  description 
of  the  amount  of  each  gas  contained  in  the  mixture  from  percent  to  fraction  and  stores  in  the  Ns() 
array  the  number  density  of  each  species.  The  order  of  storage  in  the  Ns()  array  is  N2,  C02,  and 
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He, 


3.1.3  C-Matrix  Elements  This  section  of  the  program  calculates  the  values  needed  to  form 
the  C-matrix  of  equation  (40).  The  actual  loading  of  this  matrix  is  dealt  with  in  the  next  section. 
The  subroutine  called  Interp  is  used  here.  An  explanation  and  example  of  its  use  is  given  in  the 
next  section. 

The  first  item  to  be  calculated  in  this  section  is  mjs  of  equation  (39).  This  represents  the  value 
of  the  nearest  integer  to  C;«/Ae  where  (j,  is  the  threshold  energy  loss  in  electron  volts  associated 
with  the  j***  internal  energy  state  of  species  s  and  At  is  the  bin  width  in  electron  volts.  This 
information  is  stored  in  the  array  mjs%(J,I),  where  the  symbol  %  at  the  end  of  the  array  name 
indicates  that  the  values  stored  in  this  array  are  integers.  The  index  J  refers  to  species  ( J  =  1  means 
N2,  J  =  2  means  C02)  and  I  refers  to  a  particular  internal  energy  state  (for  example  a  particular 
vibrational  state  of  N2).  One  main  and  one  nested  loop  each  are  used  for  this  calculation.  The  main 
loop  steps  through  the  individual  species  and  the  nested  loop  steps  through  the  internal  energy 
states  of  each  species.  The  function  CINT  is  internal  to  QuickBasic  and  provides  as  its  output  the 
integer  nearest  its  argument. 

The  second  part  of  this  section  of  Bolts  calculates  the  rates  used  in  the  C-matrix.  Two  main 
loops  are  used  in  this  section.  The  first  main  loop  is  for  the  calculation  of  inelastic  and  superelastic 
rates.  The  second  main  loop  is  for  the  elastic  rates.  Within  the  first  main  loop,  which  uses  the 
index  K  to  step  through  the  individual  species  of  gas,  are  two  levels  of  nested  loops.  The  first  level 
of  nesting  uses  the  index  I  to  step  through  each  of  the  individual  internal  energy  states  of  each  gas. 
This  loop  is  controlled  by  the  information  stored  in  the  jin()  array.  The  second  level  of  nesting 
contains  two  loops.  The  first  loop  uses  the  temporary  storage  arrays  TempeV()  and  TepmQj()  to 
store  energy  and  inelastic  cross  section  data  for  a  particular  internal  energy  state  of  a  particular 
gas  contained  in  the  master  three  dimensional  arrays  eV()  and  QNj().  The  second  of  these  loops 
then  uses  these  temporary  arrays  as  arguments  of  the  subprogram  Interp  in  order  to  get  a  good 
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estimate  of  the  inelastic  cross  section  at  a  particular  energy  (the  right  hand  edge  of  the  energy  bin 
under  consideration).  The  second  loop  steps  through  each  of  the  energy  bins  as  it  computes  and 
stores  the  rates.  The  array  Rjs()  is  used  to  store  the  inelastic  rates  which  have  been  calculated  for 
a  certain  species  in  a  particular  internal  energy  state  for  each  energy  bin.  The  array  RTjs()  is  used 
to  store  the  sum  of  the  individual  rates  over  all  the  internal  states  of  all  the  species  for  a  particular 
energy  bin.  This  loop  does  likewise  for  the  superelastic  rates  which  it  stores  in  the  arrays  RSjs() 
for  individual  and  RTSjs()  for  total  rates.  The  actual  calculation  of  the  inelastic  and  superelastic 
rates  are  as  shown  in  equation  (39).  The  purpose  of  the  two  IF  . .  .THEN  statements  within  this 
loop  is  to  allow  the  superelastic  rates  for  C02(  where  K=2)  to  be  calculated  with  the  user  supplied 
vibrational  temperatures  for  C02  modes  (100),  (010)  and  (001). 

The  last  part  of  this  section  concerns  calculation  of  the  elastic  rates.  The  calculations  are  as 
shown  in  equation  (39).  This  section  consists  of  just  one  FOR  . .  .NEXT  loop  which  steps  through 
each  energy  bin  with  the  index  I.  The  loop  is  controlled  by  the  variable  kc%  (number  of  bins)  which 
is  supplied  by  the  user  through  the  main  program  C020SC.  The  subprogram  Interp  is  called  three 
times  each  loop  in  order  to  provide  an  interpolated  value  of  the  momentum  transfer  cross  section 
for  each  of  the  three  gases  at  each  energy  bin.  The  variables  nup  and  nub  represent  /N  and 
respectively  from  equation  (39). 

S.I.4  (I-hC)  Matrix  This  section  of  the  program  simultaneously  loads  the  C-matrix,  mul¬ 
tiplies  it  by  the  time  step  h,  and  then  subtracts  this  resulting  product  from  the  identity  matrix. 
Combining  these  steps  allows  for  a  reduction  in  the  number  of  loops  which  would  otherwise  be 
required.  The  result  of  this  is  a  shorter  run  time  for  the  program.  The  first  part  of  this  section 
loads  into  the  a()  and  b()  arrays  the  boundary  conditions  as  specified  in  section  II. 3.  The  C-matrix 
is  then  used  to  store  the  values  of  (1-hC).  Most  of  the  values  can  be  loaded  through  a  series  of  loops. 
Those  which  must  be  loaded  outside  of  loops  are  calculated  and  loaded  first.  All  calculations  in 
this  section  are  in  accordance  with  equations  (39)  -  (40). 
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The  second  part  of  this  section  consists  of  a  single  loop  which  loads  into  the  (I-hC)  matrix 
contributions  along  the  three  diagonals  which  are  above,  along  and  below  the  main  diagonal.  These 
contributions  come  from  the  a(),  b(),  Rjs()  and  PTjs()  arrays.  Up  through  this  point  in  the  program 
Boltz,  all  contributions  will  have  been  along  one  of  these  three  diagonals. 

The  third  part  of  this  section  allows  for  contributions  to  bands  other  than  the  three  mentioned 
above.  This  part  has  one  main  loop  and  two  levels  of  nested  loops,  one  loop  in  each  level.  These 
three  levels  of  loops  step  through  each  species  of  gas  (N2  and  then  C02,  there  is  no  inelastic 
contribution  from  He),  each  energy  bin,  and  e2M:h  internal  state  of  each  species.  In  this  manner,  all 
of  the  contributions  to  the  (I-hC)  matrix  from  the  three  dimensional  Rjs()  and  RSjs()  matrices  are 
accounted  for,  in  accordance  with  equations  (39)  and  (40).  Two  IF. .  .THEN  statements  are  used 
in  the  last  nested  loop.  These  statements  provide  for  implementation  of  the  boundary  conditions 
as  applied  to  the  inelastic  and  superelastic  contributions  to  the  (I-hC)  matrix.  When  this  section 
of  Boltz  is  complete,  a  message  is  printed  to  the  screen  so  as  to  inform  the  user  of  the  progress 
Boltz  has  made  thus  far. 

3.1.5  Inversion  of  (I-hC)  Boltz  now  calculates  the  inverse  of  the  (I-hC)  matrix  (which  has 
been  stored  in  the  c()  array  in  Boltz)  in  two  parts.  The  first  performs  the  L-U  decomposition  and  the 
second  part  the  backsubstitution  described  in  Appendix  E.  Both  routines  were  taken  from  reference 
10,  where  they  were  published  in  Fortran  as  subroutines.  They  were  rewritten  in  QuickBasic  for 
use  as  a  part  of  Boltz.  They  are  used  as  part  of  the  main  code  in  Boltz,  not  as  subroutines.  When 
the  inverse  of  the  (I-hC)  matrix  has  been  completed,  the  inverse  of  (1-hC)  is  stored  in  the  c()  array 
and  a  message  is  printed  to  the  screen  informing  the  user  of  the  progress  of  Boltz. 

3.1.6  liemiive  Calenlaiion  of  Electron  Nnmber  Density  With  the  inversion  of  (I-hC)  com¬ 
pleted,  Boltz  is  ready  to  iteratively  calculate  the  normalized  electron  number  density  using  equation 
(42).  Before  the  first  iteration  can  be  calculated,  an  initial  first  guess  is  needed.  Bolts  supplies  this 
initial  guess  in  the  form  of  a  straight  line  (constant  value)  distribution.  All  elements  of  the  NO() 
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array  are  set  equal  to  a  value  of  1  before  control  is  passed  to  the  main  iteration  loop.  The  iteration 
part  of  Boltz  consists  of  a  main  loop  and  two  levels  of  nested  loops.  The  main  loop  controls  the 
maucimum  number  of  iterations  that  Boltz  may  perform.  At  present,  this  is  set  to  300.  Of  course, 
should  the  convergence  criteria  be  met  before  the  iteration  number  300,  control  will  pass  from  the 
iteration  loop  to  the  next  part  of  Boltz.  The  first  nested  loop  of  this  section  also  contains  the 
second  nested  loop.  These  two  loops  perform  the  multiplication  of  the  inverted  (I-hC)  matrix  by 
the  last  iteration  of  the  NO()  vector.  The  vector  result  of  this  operation  is  stored  in  the  NOl() 
vector.  The  latest  iteration  of  the  un-normalized  electron  number  density  is  thus  contained  in 
NOl()  at  this  point.  It  must  be  emphasized  that  what  has  actually  been  calculated  by  Boltz  is  the 
total  electron  number  density  within  each  bin,  and  not  for  a  particular  energy.  Before  the  vectors 
NO()  and  NOl()  can  be  compared  to  determine  if  the  convergence  criteria  has  been  met,  NO()  and 
N01()  must  both  be  normzdized  such  that  the  sum  of  the  elements  of  each  is  equal  to  1.  This  will 
fulfill  the  normalization  condition  of  equation  (14).  The  variables  “sum”  and  “suml”  are  used  to 
store  the  values  needed  to  normalize  NO()  and  N01()  respectively.  These  values  are  determined 
by  summing  all  of  the  values  contained  in  NO()  (for  sum)  and  N01()  (for  suml).  Each  element 
of  these  two  arrays  is  then  divided  by  its  respective  normalizing  factor.  The  values  of  NO()  and 
NOl()  are  then  compared  for  each  energy  bin  and  the  maximum  relative  difference  stored  in  the 
variable  cmax.  This  value  is  then  compared  to  the  convergence  criteria  stored  in  the  variable  conv. 
If  cmax  is  less  than  conv,  then  the  main  iteration  loop  is  terminated.  If  the  convergence  test  fails, 
the  NO()  vector  is  updated  with  the  values  from  NOl()  and  the  iteration  procedure  is  repeated. 

Once  the  desired  level  of  convergence  is  achieved,  NO()  is  updated  with  the  last  normalized 
values  from  NO(),  and  plotting  occurs.  Since  NO()  is  a  double  precision  array  and  the  arguments 
of  Plotl  are  single  precision,  NO()  must  be  temporarily  stored  in  the  single  precision  array  PNO() 
for  plotting. 


S.1.7  Transport  Coefficients  Once  the  normalized  electron  number  density  hew  been  deter¬ 
mined,  Boltz  calculates  the  transport  coefficients  required  for  output  to  C020SC.  The  first  step  in 
this  section  of  the  program  is  to  calculate  the  values  of  St  and  5*  required  in  equation  (56).  These 
values  are  stored  in  the  abar()  and  bbar()  arrays. 

The  first  transport  coefficient  calculated  is  the  drift  velocity.  Boltz  accomplishes  this  using 
equation  (59)  as  the  basis.  Because  the  drift  velocity  calculation  is  performed  in  units  of  S.I.,  the 
value  which  Boltz  reports  is  divided  by  .01  so  as  to  convert  the  output  to  cm/s.  Thus  the  reported 
value  is  in  units  most  often  found  in  the  literature.  Once  the  drift  velocity  has  been  computed, 
Boltz  reports  its  value  to  the  screen  along  with  many  of  the  original  parameters  that  went  into  the 
calculation.  This  enables  the  user  to  press  SHIFT-PRINT  SCREEN  and  have  all  this  information 
sent  to  the  printer  at  once. 

Boltz  next  calculates  the  inelastic  excitation  rates  for  the  eight  vibrational  levels  of  N2  and 
modes  (100)  and  (001)  of  C02.  These  calculations  are  accomplished  using  equation  (61).  For  N2, 
the  inelastic  excitation  rates  for  each  of  the  eight  individual  levels  and  their  total  are  reported  to 
the  screen.  For  C02,  the  inelastic  excitation  rates  for  (100)  and  (001)  modes  are  reported.  Boltz 
also  calculates  a  weighted  total  for  N2  and  stores  this  in  TOTN2R.  This  weighted  sum  is  used  in 
the  calculation  of  the  pumping  term  for  N2.  An  approximation  is  made  here  by  assuming  that  the 
eight  vibrational  levels  of  N2  are  equally  spaced.  This  means  that  an  N2  molecule  excited  to  the 
v=2  state  can  vibrationally  excite  two  C02  molecules,  in  the  v=3  state  an  N2  molecule  can  excite 
3  C02  molecules,  etc. 


S./.S  Pumping  Terms  The  objective  of  Bolts  is  to  supply  C020SC  with  the  pumping  terms 
Wa,  Wb  and  Wc.  The  formula  as  required  by  C020SC  is 


W 

derit  ■  (1/fcav) 


=  Pump  Rate  (sec“  * ) 


(50) 
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where  W  can  be  Wa,  Wb  or  Wc  and  (3:162) 

Wa  =  N fNcoiXcoiiooi) 

Wb  =  NtNco2Xcot(ioo) 

Wc  =  NtNN‘iXS2 

where  Ng  is  the  electron  number  density  (unnormalized)  and  x  is  the  respective  inelastic  excitation 
rate.  To  calculate  Ng,  Boltz  uses  (3:167) 

Ng  =  (current  density)/(drift  velocity  x  electron  charge) 

This  allows  Boltz  to  complete  the  calculation  of  parameters  to  be  passed  b2tck  to  C020SC.  The 
final  calculation  performed  by  Boltz  is  energy  balemce.  This  acts  as  a  check  on  the  quality  of  the 
calculations  performed  to  arrive  at  the  electron  number  density. 

3.2  Subprogram  Inierp 

Interp  is  a  subprogram  used  by  Boltz  to  perform  linear  interpolation.  This  is  important 
because  the  cross  section  data  files  contain  data  only  at  certain  values  of  energy.  An  assumption  is 
made  here  that  linear  interpolation  will  provide  a  sufficiently  accurate  estimate  for  a  cross  section 
value  which  falls  somewhere  between  two  known  data  points.  This  is  based  upon  the  assumption 
of  linear  behavior  of  the  cross  sections  between  these  known  points.  Before  Interp  was  written  a 
cubic  spline  based  interpolation  program  was  tried.  Due  to  the  very  abrupt  changes  in  the  cross 
section  data  curves  (14:62-67),  the  cubic  spline  results  tended  to  overshoot  and  give  interpolated 
cross  section  values  far  removed  from  the  known  data.  In  some  instances  negative  values  for  the 
cross  section  estimates  were  output  by  the  cubic  spline  routine.  Therefore,  a  linear  interpolation 
based  routine  was  written.  A  copy  of  the  listing  for  subprogram  Interp  is  in  Appendix  C.  Interp  is 
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based  in  part  on  the  subroutine  SPLINT  in  reference  10. 


There  are  five  parameters  in  the  argument  of  Interp.  The  first  four  are  input  parameters 
supplied  to  Interp  by  the  calling  program,  and  the  last  is  the  output  Interp  provides  to  the  calling 

program.  These  parameters  are; 

x()  -  an  array  with  the  known  values  for  the  abscissas 

y()  -  an  array  with  the  known  ordinates  corresponding  to  the  x()  array 

N  -  total  number  of  points  contained  in  the  x()  array 

xin  -  the  input  abscissa  value  for  which  an  ordinate  value  is  to  be  determined  by  Interp 
yout  -  the  output  interpolated  value 


The  first  action  Interp  performs  is  to  see  if  the  input  abscissa  value  is  within  the  range  for 
which  non-zero  cross  section  values  exist.  If  it  is  outside  this  range,  this  is  sensed  by  the  IF. . .  THEN 
and  ELSEIF  statements.  These  statements  direct  Interp  to  set  yout  equal  to  zero  and  return  to 
the  main  program. 


If  Interp  has  determined  that  a  non-zero  cross  section  value  exists  for  xin,  it  then  proceeds  to 
find  one.  Interp  uses  the  vririables  klo  and  khi  as  place  markers  for  the  low  and  high  values  in  the 
array  x().  Interp  initializes  these  place  markers  with  the  lowest  and  highest  index  values  of  the  x() 
array  (1  and  N).  Using  an  IF. .  .THEN. . .  ELSE  construction,  Interp  performs  a  looping  bisection 
search.  Interp  then  denotes  the  lower  value  as  x(klo)  and  the  higher  value  as  x(khi). 


Now  that  Interp  knows  the  two  abscissas  between  which  xin  lies,  it  is  ready  to  try  to  obtain 
an  ordinate  to  assign  to  xin.  First  though,  Interp  performs  an  internal  check  to  make  sure  that  it 
has  not  accidentally  assigned  the  same  point  in  the  x()  array  to  both  x(klo)  and  x(khi).  Once  the 
internal  error  check  is  complete,  Boltz  determines  the  value  of  yout  to  assign  to  xin  based  upon 


the  slope  between  the  points  (x(klo),  y(klo))  and  (x(khi),  y(khi)).  Boltz  assigns  to  the  variable  h 
the  difference  in  abscissa  values  for  these  two  points.  Even  though  the  variable  h  is  used  elsewhere 
in  Boltz,  this  does  not  cause  a  problem.  Variables  defined  within  a  subprogram  are  local  to  that 
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subprogram.  This  is  an  advantage  of  defining  Interp  as  a  subprogram,  instead  of  as  a  subroutine 
(in  which  variables  become  global).  Boltz  then  uses  h  in  the  formula 

yout  =  y(klo)  +  ((xin  -  x(klo))/h)  ♦  (y(khi)  -  y(klo)) 

The  interpolated  ordinate  value  Interp  has  assigned  to  xin  is  then  passed  back  to  Boltz  as  the 
parameter  yout. 

3.3  Plotting  Routines 

The  two  plotting  routines  used  are  Plotl  and  Plot2.  The  listings  for  these  subprograms  are 
in  Appendix  C.  Plotl  is  used  by  Boltz  to  plot  the  electron  number  density  and  Plot2  is  used  by 
C020SC  to  plot  laser  output  power  and  output  total  energy.  Both  subprograms  provide  linear  and 
semilog  plots.  Plot2  has  the  additional  capability  of  allowing  the  user  to  specify  the  minimum  and 
maximum  values  for  both  coordinate  axes  and  to  replot  a  particular  power  or  energy  curve  as  many 
times  as  desired.  Both  routines  make  use  of  the  subprogram  cgadump  (16:156),  which  can  speed 
up  printing  out  a  plot  on  a  printer  (as  opposed  to  using  SHIFT-PRINT  SCREEN).  Since  Plotl 
and  Plot2  are  nearly  identical,  their  descriptions  are  combined  here,  with  differences  highlighted  as 
they  occur. 

3.3.1  Input  Parameters  Ail  of  the  parameters  in  the  plotting  subprograms  are 

supplied  as  input  by  the  calling  program.  For  Plotl  these  are: 
de  -  abscissa  coordinate  spacing 
kc%  -  number  of  points  to  be  plotted 
NO()  -  ordinate  value  to  be  plotted 

h,  T,  P,  EN  -  numerical  values  of  time  step,  ambient  temperature,  pressure,  and  E/N  to  be 
displayed  as  information  on  the  plots 
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For  Plot2  the  parameters  are: 

hi,  h2  -  integration  time  steps  (abscissa  coordinate  spacing  in  units  of  tcav) 
jChStep  -  the  point  in  the  ordinate  data  file  (the  index  number  of  that  array  point)  where 
integration  time  steps  are  changed 

ChTimeStep  -  if  equal  to  zero  then  integration  time  steps  are  not  changed  (only  hi  is  used), 
and  if  equal  to  other  than  zero  then  time  steps  change  from  hi  to  h2  at  jChStep 

tCav  -  unit  of  cavity  lifetime  as  determined  by  C020SC  and  used  in  Plot2  to  provide  actual 
time  values  for  the  abscissa  coordinates 

jmax  -  number  of  points  in  the  ordinate  array  to  be  plotted 

Yquant()  -  ordinate  array  to  be  plotted 

titles  -  character  string  to  be  displayed  on  the  plot 

3.3.S  Linear  Plot  Lineal  and  semilog  plots  are  generated  using  these  routines.  The  code 
used  to  produce  each  type  is  very  nearly  the  same.  The  linear  plot  code  will  be  discussed  first.  The 
semilog  discussion  will  only  highlight  differences  between  the  two. 

3.3.3  Axes  and  Tick  Marks  The  first  step  of  the  actual  plotting  is  to  draw  the  horizontal 
and  vertical  axes  and  their  corresponding  tick  marks.  Before  these  lines  are  drawn  however,  the 
plotting  routines  first  determine  numerical  ranges  over  which  they  will  be  working  for  both  axes. 
In  the  case  of  Plotl,  the  minimum  abscissa  coordinate  is  the  right  hand  edge  of  the  first  bin  and 
the  maximum  is  the  right  hand  edge  of  the  last  bin.  These  values  are  stored  in  XMin  and  XMax 
respectively.  A  FOR. . .  NEXT  loop  is  then  performed  so  as  to  sift  through  the  NO()  array  and  store 
its  minimum  and  maximum  values  in  the  variables  Min  and  Max.  Plot2  can  also  automatically 
determine  the  range  of  values  to  be  plotted.  For  the  first  Plot  during  a  particular  call,  Plot2  will 
branch  to  the  label  “autocoords” .  The  minimum  and  maximum  values  are  determined  in  a  similar 
manner  to  that  of  Plotl.  The  difference  is  that  the  minimum  abscissa  coordinate  is  zero  and  the 
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maximum  corresponds  to  the  time  for  the  last  element  of  Yquant().  If  there  is  no  change  in  time 
step,  the  maximum  will  be 

XMax  =  (hi  •  tcav  •  10"®)  ■  jmax 

This  also  converts  the  time  associated  with  jmax  from  units  of  tcav  to  nanoseconds.  If  a  change  in 
time  step  does  occur,  then  the  maximum  will  be 

XMax  =  (jChStep  hi  +  (jmax  —  jChStep)  •  h2)  tcav  ■  10"® 

Later  in  Plot2,  the  user  may  choose  to  replot  a  particular  set  of  data  and  is  then  given  the 
choice  of  manually  setting  these  minimum  and  maximum  values.  This  occurs  in  the  section  with  the 
label  “mancoords” .  After  inputting  the  range  from  the  abscissa  coordinates,  the  user  is  offered  the 
choice  of  letting  Plot2  automatically  set  the  range  for  the  ordinate  values  or  enter  them  manually. 
Manual  entry  affords  the  user  the  opportunity  to  construct  plots  from  different  computer  runs 
which  have  their  data  displayed  over  the  same  ranges  for  ease  of  comparison. 

In  preparation  of  drawing  on  the  screen,  the  routines  first  clear  the  screen  with  CLS  0  com¬ 
mand.  The  resolution  is  then  set  to  640x200  with  the  SCREEN2  command.  The  WINDOW  and 
VIEW  commands  ensure  that  any  previous  graphics  commands  which  might  affect  the  display  are 
negated.  The  axes  and  tick  marks  are  then  drawn  using  the  LINE  command  and  two  FOR. . .  NEXT 
loops. 


3.3.4  Label  Tick  Marks  In  this  section  of  the  code,  tick  marks  have  numerical  values  printed 
next  to  them  and  labels  are  printed  in  the  screen.  The  label  identifying  the  plot  is  printed  first 
along  with  information  concerning  input  parameters  which  may  have  affected  the  plot.  The  range 
of  abscissa  and  ordinate  values  is  calculated  and  stored  in  the  variables  XRange  and  YRange.  The 
increment  between  abscissa  and  ordinate  tick  marks  is  also  calculated  and  stored  in  XIncr  and 
YIncr.  These  incremental  values  are  then  used  in  the  FOR. .  .NEXT  loops  which  label  each  tick 
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mark. 


3.3.5  Plot  the  Data  Before  the  routines  draw  lines  connecting  data  points,  the  pixel  loca¬ 
tions  on  the  screen  must  be  redefined  in  order  to  accommodate  the  axes.  If  this  is  not  done,  then 
the  axes,  tick  marks  and  plot  will  bear  no  relationship  to  each  other.  The  area  where  the  plotting 
will  actually  occur  is  defined  with  the  VIEW  command,  the  arguments  of  which  define  the  pixel 
coordinate  location  of  the  lower  left  and  upper  right  hand  corners  of  the  plot.  This  lower  left  hand 
corner  is  set  to  the  point  where  the  two  axes  meet.  The  upper  right  hand  corner  is  set  to  the 
corner  of  the  screen.  The  pixels  themselves  are  scaled  using  the  WINDOW  command.  This  allows 
pixel  coordinates  to  be  specified  according  to  an  arbitrary  scale.  This  means  that  during  plotting, 
QuickBasic  will  automatically  convert  the  values  of  the  coordinates  of  the  points  to  be  plotted  to 
pixel  coordinates  and  then  map  that  point  onto  the  screen.  This  eliminates  the  need  to  set  up  a 
loop  to  accomplish  these  tasks  manually. 

With  the  viewing  port  and  pixel  scading  accomplished,  the  routines  can  begin  plotting.  In 
Plotl  this  is  accomplished  with  a  single  FOR.. .NEXT  loop.  The  abscissa  and  ordinate  values  of  the 
points  are  sequentially  stored  in  the  four  variables  (xl,yl)  and  (x2,y2)  and  a  line  drawn  between 
them.  In  PIot2,  if  no  change  of  integration  time  step  has  occurred,  then  the  plotting  procedure 
is  the  same  as  that  in  Plotl.  When  a  change  of  integration  time  step  has  been  used,  control  in 
Plot2  is  then  transferred  to  the  point  with  the  label  “varplot” .  Plotting  then  occurs  as  a  two  step 
process.  Points  which  occur  in  time  before  the  change  are  plotted  with  a  spacing  based  upon  hi. 
Those  which  occur  after  the  change  are  plotted  with  a  spacing  based  upon  h2.  For  the  case  where 
one  point  is  before  the  change  and  the  other  after,  Plot2  applies  the  correct  time  factor  to  each 
and  plots  them. 

When  linear  plotting  has  finished,  the  user  is  offered  the  opportunity  to  send  the  plot  to  the 
printer.  This  dump  of  the  plot  is  handled  by  the  subprogram  cgadump.  This  routine  was  taken 
from  Cooper  (16:156).  On  an  Apple  IIGS  with  an  Applied  Engineering  PC  Transporter  installed. 
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using  SHIFT-PRINT  SCREEN  to  dump  a  plot  to  the  printer  required  4  minutes  45  seconds  to 
print  the  plot.  Using  cgadump  reduced  this  to  2  minutes  30  seconds. 

The  user  of  Plot2  is  then  given  the  chance  to  replot  this  set  of  data.  Choosing  to  do  this  will 
cause  control  to  branch  back  to  the  label  “stplot”,  and  the  minimum  and  maximum  values  of  the 
coordinate  axes  may  then  be  set  manually. 

3.3.6  Log  Plot  Upon  completion  of  the  linear  plot,  the  routines  perform  a  plot  of  the  same 
data,  only  this  time  with  the  ordinate  axis  in  common  log  scale.  This  part  of  the  routines  is  nearly 
identical  with  the  linear  plot  section  of  the  code.  The  first  change  is  that  the  ordinate  values  must  be 
converted  to  common  log  values.  Since  QuickBasic’s  internal  LOG  command  performs  computation 
of  a  natural  logarithm  only,  these  log  value  are  multiplied  by  .434294482.  This  effects  conversion 
from  natural  to  common  logarithm.  The  second  change  occurs  in  determining  the  minimum  and 
maximum  abscissa  values.  For  ease  of  interpretation  of  the  log  scale  used  on  the  ordinate  axis, 
the  minimum  and  maximum  values  are  converted  to  integer  values.  This  is  accomplished  using  the 
INT  command. 

3.4  Data  Files 

The  data  files  are  saved  on  disk  under  the  names  N2DATA,  C02DATA  and  HEDATA.  These 
contain  the  cross  section  data  for  N2,  C02  and  He.  All  are  saved  as  sequential  ASCII  files.  This 
method  of  storage  was  chosen  since  it  would  provide  easy  access  to  the  data  files  for  most  users, 
should  the  data  file  need  to  be  modified.  Any  word  processor  capable  of  reading  an  ASCII  file  may 
be  used  to  change  the  data  in  these  three  files.  It  is  not  necessary  that  the  user  have  a  copy  of 
Microsoft  QuickBasic  to  do  this.  The  three  data  files  are  in  Appendix  D. 

The  N2  and  C02  data  files  are  essentially  the  same.  The  first  record  contains  one  field.  This 
field  is  stored  in  Boltz  in  the  array  jin().  It  tells  how  many  internal  energy  states  are  contained  in 
the  data  file.  The  second  record  contains  three  fields.  The  first  field  of  this  record  tells  how  many 
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records  there  are  for  this  particular  internal  energy  state.  The  second  field  contains  the  name  of  the 
gas  and  the  internal  energy  state.  The  third  field  is  the  threshold  inelastic  energy  loss  for  this  state, 
and  is  the  best  means  for  positively  identifying  to  as  which  group  of  data  in  reference  14  this  data 
set  belongs.  The  rest  of  the  records  for  that  internal  energy  state  will  consist  of  two  fields  each. 
The  first  field  will  contain  an  energy  in  electron  volts  and  the  second  its  corresponding  inelastic 
cross  section  in  cm^,  from  reference  14.  After  the  last  record  for  a  particular  internal  energy  state, 
will  be  the  first  record  of  the  next  internal  state.  This  pattern  will  continue  until  all  of  the  inelastic 
cross  section  data  is  exhausted,  and  the  momentum  transfer  cross  section  data  begins. 

.\11  three  files  share  the  same  basic  set  up  for  the  momentum  transfer  cross  section  data.  For 
N2  and  C02,  this  follows  immediately  after  the  inelastic  record.  For  He,  this  is  the  only  data  in 
the  file.  The  first  record  consists  of  two  fields.  The  first  field  tells  how  many  momentum  transfer 
cross  section  records  are  contained  in  the  file.  The  second  field  is  the  name  of  the  gas.  The  rest  of 
the  records  consist  of  two  fields  each.  The  firs  field  is  the  energy  in  electron  volts  and  the  second 
its  corresponding  momentum  transfer  cross  section  in  cm*. 

3.5  Interface  wtth  C020SC 

Boltz  and  its  accompanying  subprograms  were  written  so  as  to  require  only  minimal  changes 
to  C020SC.  Users  wishing  to  modify  C020SC  need  to  be  aware  however,  of  how  Boltz  fits  into 
C020SC  and  its  special  requirements.  Failure  to  allow  for  the  requirements  of  Boltz  could  lead  to 
a  fat2tl  error  in  attempting  to  run  C020SC  or  cause  erroneous  output  data  to  be  generated.  Listed 
below  are  all  the  places  in  C020SC  where  an  impact  on  the  main  level  code  from  Boltz  can  be 
seen.  The  listings  for  these  portions  of  C020SC  are  in  Appendix  B. 

One  of  the  most  important  considerations  for  using  Boltz  with  C020SC  is  memory  require¬ 
ments.  Depending  on  the  number  of  bins  specified  by  the  user,  Bolts  could  temporarily  need  all  of 
the  remaining  memory  in  the  computer.  A  640K  computer  can  handle  150  bins  when  running  the 
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compiled  version  of  C020SC,  and  93  bins  when  run  from  the  QuickBasic  shell.  Because  Boltz  can 
be  very  memory  intensive,  it  should  be  run,  if  at  all  possible,  before  any  other  subprogram.  Before 
Boltz  terminates  and  returns  control  to  C020SC,  it  de-allocates  all  unnecessary  storage. 

At  the  beginning  of  C020SC  are  the  declaration  statements  for  the  subprograms.  Five  of 
these  subprograms  pertain  to  Boltz  or  were  added  to  C020SC  as  part  of  this  project.  Boltz  and 
Plot2  are  subprograms  which  aire  called  from  the  main  level  code  of  C020SC.  Interp  and  Plotl 
are  called  from  Boltz  and  cgadump  is  called  from  Plotl  and  Plot2.  After  the  last  DECLARE 
statement  is  a  CLEAR  statement.  This  was  added  to  C020SC  so  as  to  allocate  more  space  to 
the  area  of  memory  called  the  stack.  This  was  done  to  accommodate  the  extra  space  needed  in 
the  stack  by  the  newly  added  subprograms.  If  more  subprograms  are  added,  the  memory  block  set 
aside  by  the  CLEAR  command  may  have  to  increase.  Unfortunately,  stack  space  competes  with 
memory  requirements  for  string  space.  Setting  aside  too  much  memory  for  the  stack  may  deprive 
some  other  portion  of  the  program  memory,  such  as  string  variables,  space  it  needs  to  run.  One 
possible  remedy  to  this  problem  is  to  run  C020SC  as  a  compiled  program,  rather  than  from  the 
QuickBasic  shell.  This  will  free  up  the  memory  normally  occupied  by  the  QuickBasic  shell  for  use 
by  C020SC,  which  is  why  the  compiled  version  can  be  run  with  a  higher  bin  setting  than  that 
available  when  C020SC  is  run  from  the  QuickBasic  shell. 

The  nominal  input  parameter  list  has  been  changed  to  include: 

EN  -  E/N  in  volts  cm2 

VT()  -  vibrational  temperature  for  C02  modes  (100),  (010)  and  (001) 
kc%  -  number  of  energy  bins 
de  -  bin  width  in  electron  volts 
jd  -  plasma  tube  current  density  in  amps/cm2 

Prev$  -  a  string  variable  used  to  tell  C020SC  if  the  user  has  decided  to  use  the 
previously  computed  kinetics  (thereby  negating  the  need  to  call  Boltz  again) 
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These  same  variables  are  again  used  in  the  main  level  code  of  C020SC  in  the  sections  labeled 
“inputroutine”  and  SELECT  CASE  ChoiceS.  These  two  sections  of  the  code  are  used  to  display 
to  the  user  the  current  default  settings  used  by  C020SC  and  to  allow  the  user  the  opportunity  to 
change  these  settings. 

In  the  section  labeled  “pump  rates  and  effective  spontaneous  emission  rate”  is  where  Boltz 
is  actually  called.  An  IF. .  .THEN  statement  is  used  to  check  if  the  user  has  decided  to  use  Boltz 
or  the  kinetic  data  from  the  previous  run.  After  Boltz  is  run,  a  check  is  made  to  see  if  the  user 
is  satisfied  with  the  Boltzmann  calculation  and  wishes  to  continue  the  program  or  return  to  the 
main  menu  and  change  one  or  more  parameters  (for  example  bin  width).  If  the  user  has  chosen  to 
continue  the  program,  a  message  indicating  that  the  laser  portion  of  the  calculation  is  in  progress 
is  sent  to  the  screen.  C020SC  then  uses  the  parameters  passed  to  it  from  Boltz  to  calculate  the 
pump  rates  to  be  used  by  C020SC  in  its  rate  equations. 

The  final  change  to  C020SC  is  in  the  section  labeled  “Output  routines”.  Here,  Plot2  has 
been  added.  Each  call  to  Plot2  is  preceded  by  storage  in  the  string  variable  titleS.  This  is  a  title 
to  be  displayed  on  the  plot  produced  by  Plot2. 

S.6  User’s  Guide 

Boltz  offers  the  user  of  C020SC  the  ability  to  account  for  the  effects  of  the  electric  discharge 
excitation  in  computing  C02  laser  output.  However,  it  also  places  new  demands  on  the  user. 
Specifically,  the  user  must  know  how  to  activate  and  use  the  new  features  of  C020SC  provided  by 
Boltz.  This  section  will  cover  the  new  features  of  C020SC  and  explain  how  to  use  them. 

The  procedure  for  booting  C020SC  is  as  before.  However,  the  first  screen  presented  to  the 
user  has  eight  new  input  parameters.  These  new  input  parameters  are  labeled  y  through  y7.  The 
results  section  explores  the  effects  of  these  various  parameters  and  explains  the  reasons  why  the 
nominal  values  provided  in  C020SC  were  chosen. 
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The  first  parameter,  y,  allows  the  user  to  select  a  value  for  E/N.  This  is  the  ratio  of  the 
electric  field  strength  in  volts/cm  to  the  total  neutral  gas  number  density  in  cm~^.  Therefore  the 
units  for  E/N  are  Volts  cm^. 

The  next  three  parameters,  yl  through  y3,  allow  the  user  to  specify  the  vibrational  tempera^ 
tures  of  C02  modes  (100),  (010),  and  (001).  The  unit  for  these  vibrational  temperatures  is  degrees 
kelvin. 

The  fifth  parameter,  y4,  is  the  number  of  energy  bins.  This  is  the  parameter  that  decides 
into  how  many  cells  or  bins  the  energy  axis  will  be  divided.  The  more  bins  used,  the  better  the 
accuracy  will  usually  be  for  a  particular  run.  There  is  an  upward  limit  of  150  bins  if  C020SC  is 
run  from  a  compiled  version,  or  93  bins  if  run  from  the  QuickBasic  shell.  The  user  must  be  aware 
that  a  larger  number  of  bins  can  dramatically  affect  the  run  time  for  C020SC.  As  a  guide,  an 
8MHz  computer  with  a  math  coprocessor  will  perform  a  20  bin  calculation  in  about  30  seconds, 
and  a  100  bin  calculation  in  about  10  minutes. 

The  sixth  parameter,  y5,  is  the  width  in  units  of  electron  volts  of  an  individual  energy  bin. 
Each  bin  has  the  same  width.  Multiplying  the  number  of  bins  by  the  bin  width  gives  the  total 
energy  range  over  which  the  Bolt*  calculations  are  performed. 

The  seventh  parameter,  y6,  is  the  current  density,  expressed  in  units  of  amps/cm*.  This 
is  actually  a  circuit  parameter  provided  by  the  user.  It  is  used  by  Bolt*,  after  the  electron  drift 
velocity  has  been  calculated,  to  determine  the  electron  number  density. 

The  eighth  parameter,  y7,  offers  the  user  the  opportunity  to  use  the  previously  computed 
electron  kinetics,  and  avoid  running  Bolt*  again.  This  will  save  time  if  the  user  has  decided  to 
perform  all  C020SC  runs  using  one  particular  set  of  variables  which  affect  the  electron  kinetics, 
and  just  vary  those  parameters  which  do  not  affect  the  electron  kinetics.  The  parameters  requiring 
that  Bolt*  be  run  again  if  changed  are:  d,  e,  f,  g,  k,  y,  yl,  y2,  y3,  y4,  yS,  y6  and  y7. 

Once  the  final  choice  of  parameters  has  been  made,  pressing  return  will  start  the  Boltzmann 
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portion  of  the  calculation.  Boltz  will  print  to  the  screen  messages  to  let  the  user  know  the  progress 
it  has  made.  When  calculation  of  the  normalized  electron  number  density  is  complete,  its  linear 
plot  will  be  displayed.  At  the  bottom  of  the  screen,  the  user  is  asked  if  a  copy  should  be  sent 
to  the  printer.  After  the  linear  plot,  a  semilog  version  is  displayed  on  the  screen  and  the  user  is 
again  afforded  the  opportunity  to  send  a  copy  to  the  printer.  After  the  plots  have  been  displayed, 
the  pertinent  input  parameters  and  calculated  transport  coefficients  are  displayed  on  the  screen. 
Pressing  RETURN  will  then  present  the  user  the  chance  to  return  to  the  main  menu  or  continue 
the  laser  calculation.  Choosing  to  continue  the  laser  calculation  causes  a  confirmation  message  to 
be  displayed  on  the  screen. 

The  only  change  to  C020SC  after  this  is  the  ability  to  now  plot  the  laser  output  power  and 
energy  as  a  function  of  time  from  within  C020SC.  The  first  plot  to  be  displayed  is  a  linear  plot  of 
power  versus  time.  The  plotting  routine  automatically  chooses  scales  for  the  axes  for  the  first  plot. 
After  this  plot  has  been  displayed  on  the  screen,  the  user  may  choose  to  plot  this  data  set  again.  If 
the  choice  is  made  to  replot  the  data,  the  user  will  then  be  given  the  choice  of  manually  selecting 
or  letting  the  program  select  the  range  over  which  each  axis  is  to  be  plotted.  After  the  linear  plot, 
a  semilog  plot  of  the  same  data  will  be  displayed,  and  the  same  choices  offered  to  the  user.  This 
same  cycle  will  then  be  repeated  for  the  plot  of  energy  versus  time. 
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IV.  Results 


The  results  presented  here  fall  into  two  categories.  First,  a  comparison  is  made  between 
predictions  from  Boltz  and  the  published  results  of  others.  Second,  a  study  is  made  of  the  effects 
on  laser  output  (power  and  energy)  as  various  parameters  related  to  Bolts  are  changed.  Some 
of  the  following  studies  were  made  using  the  nominal  parameters  of  C020SC,  while  others  had 
parameters  set  to  match  the  conditions  for  results  published  by  others.  At  the  end  of  this  section 
are  the  rationale  for  the  nominal  parameters  and  the  suggested  operating  limitations  for  C020SC, 
as  these  parameters  and  limitations  apply  to  the  Bolts  modification  made  to  C020SC. 

4.1  Comparison  with  Published  Results  for  a  Single  Gas  (NS) 

In  this  section,  comparisons  of  predictions  from  Bolts  are  evaluated.  Data  and  information 
are  also  presented  to  gain  an  understanding  of  what  happens  when  the  number  of  energy  bins  is 
varied.  This  is  an  attempt  to  help  the  user  of  C020SC  appreciate  the  tradeoffs  which  arise  when 
deciding  whether  to  use  a  large  or  small  number  of  bins.  To  simplify  the  situation,  drift  velocities 
and  excitation  rates  are  calculated  for  a  single  gas  (N2).  Two  sources  of  published  results  were 
available  for  comparison  with  Boltz.  An  article  containing  computer  predictions  similar  in  nature  to 
Boltz  was  published  by  Rockwood  and  Greene  (8:383,393).  Their  results  are  used  here  to  compare 
with  Boltz,  predictions  made  for  drift  velocity  and  inelastic  excitation  rates  for  N2.  Data  for  the 
experimental  determination  of  electron  drift  velocity  has  been  published  by  Crompton  and  Huxley 
(15:627-628).  Comparison  with  their  results  is  also  made  here. 

The  computer  runs  for  the  drift  velocity  data  presented  in  Table  1  were  made  at  a  temperature 
of  293  K.  The  computer  used  for  these  and  ail  subsequent  calculations  was  an  Apple  IIGS  with 
an  Applied  Engineering  PC  Transporter,  math  coprocessor  and  640K  of  ram  installed.  The  range 
of  E/N  for  the  drift  velocity  comparisons  is  from  five  to  forty  Townsends  (Td).  Rockwood  (8:383) 
suggested  this  range  so  as  to  restrict  the  calculations  to  the  region  where  vibrational  excitation 
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E/N 

Table  1. 

Nitrogen  Drift 

(lO-i^Vcm^) 

Crompton 

Rockwood  Error 

5 

1.095 

1.09 

0.5% 

10 

1.838 

1.78 

3.2% 

15 

2.48 

2.43 

2.0% 

20 

3.09 

3.06 

1.0% 

30 

4.17 

4.21 

1.0% 

40 

5.18 

5.26 

1.5% 

Velocity  Comparisons 
t;j(10®cm/s) 

CO2OSC(30)  Error  C020SC(100)  Error 


1.075  1.8%  1.070  2.3% 

1.816  1.2%  1.793  2.4% 

2.516  1.5%  2.46  0.8% 

3.18  2.9%  3.10  0.3% 

4.38  5.0%  4.26  2.2% 

5.44  5.0%  5.28  1.9% 


would  dominate.  That  restriction  is  followed  here  and  is  adopted  as  a  suggested  operating  limitation 
for  C020SC. 

In  Table  1,  the  number  in  parentheses  following  C020SC  is  the  number  of  bins  used  in  the 
calculation.  The  energy  range  was  0  to  2.5  eV  for  all  of  the  C020SC  calculations.  On  average,  the 
Rock  wood  calculations  do  the  best  over  the  entire  range.  The  100  bin  C020SC  calculation  comes 
in  close  behind.  The  30  bin  C020SC  calculations  would  have  been  excellent  except  for  the  large 
amount  of  error  in  the  higher  E/N  values.  This  table  gives  the  user  of  C020SC  an  idea  of  some  of 
the  potential  errors  and  magnitudes  of  the  Bolt*  portion  of  the  C020SC  program.  If  the  user  can 
live  with  these  errors,  the  30  second  run  time  of  the  30  bin  calculation  may  be  more  appealing  than 
the  9.5  minute  run  time  of  the  100  bin  calculation.  The  effect  of  changing  the  time  step  h  on  the 
computed  drift  velocity  was  also  examined.  It  was  found  that  for  h  equal  to  10~'^  to  10^^  seconds, 
the  drift  velocity  varied  about  .1%.  Time  steps  above  10^’  caused  the  computer  to  overflow.  Time 
steps  below  10~^^  resulted  in  too  many  iterations,  and  the  program  aborted  the  calculation  before 
the  convergence  criteria  was  met. 


In  addition  to  the  above  data,  Rockwood  also  reports  the  results  of  calculations  made  with 
his  program  for  N2  at  300  K.  He  used  a  250  point  energy  grid  with  a  bin  width  of  .01  eV.  This 
calculation  was  repeated  with  C020SC  using  a  100  point  grid  with  a  bin  width  of  .025  eV.  The 
results  listed  in  Table  2  are  for  electron  drift  velocity  and  inelastic  excitation  rates  for  the  eight 
vibrational  states  of  N2. 


The  drift  calculations  from  the  two  different  programs  are  in  very  good  agreement.  The 
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Table  2.  Comparison  with  Rock  wood  250  Point  Grid 


Rockwood 

C020SC 

Difference 

i;rf(10®cm/s) 
Excitation  Rates 

1.782 

1.793 

0.6% 

N2  v=l  (cm^/sec) 

3.493E-10 

3.490E-10 

0.1% 

2 

6.532E-11 

6.635E-11 

1.6% 

3 

2.755E-11 

2.824E-11 

2.5% 

4 

4.516E-12 

4.688E-12 

3.8% 

5 

8.762E-13 

9.259E-13 

5.7% 

6 

9.911E-15 

l.llOE-14 

12.0% 

7 

7.474Ei-16 

8.616E-16 

15.3% 

8 

0 

0 

0 

Weighted  Excitation 

5.851E-10 

5.899E-10 

0.82% 

weighted  excitation  is  calculated  by  multiplying  each  excitation  rate  by  its  vibrational  state  number, 
and  then  summing  the  resulting  products.  As  explained  in  section  V,  the  excitation  of  the  upper 
laser  level  of  C02  occurs  in  part  as  a  result  of  the  transfer  of  vibrational  energy  from  N2.  It  is 
assumed  here  that  the  eight  vibrational  states  on  N2  are  equally  spaced,  and  that  the  number 
of  C02  molecules  excited  by  a  particular  N2  molecule  is  directly  proportional  to  the  vibrational 
state  of  the  N2  molecule.  Hence,  the  comparison  of  the  weighted  excitations.  It  is  interesting  note 
that,  even  though  the  difference  between  the  rates  calculated  by  the  two  programs  grows  as  the 
vibrational  state  increases  in  energy,  the  actud  value  of  the  inelastic  excitation  rates  shrinks  very 
quickly  with  increasing  energy.  This  causes  the  weighted  excitation  rates  to  actually  be  much  closer 
than  might  otherwise  appear.  Because  the  excitation  rates  calculated  by  Boltz  play  a  major  role 
in  the  calculation  of  the  pump  rates  to  be  supplied  to  C020SC,  it  is  encouraging  that  they  are  in 
good  agreement  with  the  published  results  of  a  commercially  available  program. 

The  results  of  this  next  section  examine  how  the  range  of  the  energy  axis  used  affects  the 
results  from  Boltz  as  the  number  of  bins  is  held  constant  at  100.  The  information  presented  here  is 
intended  to  help  the  user  of  C020SC  decide  on  a  suitable  bin  width,  since  that  is  one  of  the  user 
definable  parameters  of  C020SC.  This  number  of  bins  was  chosen  since  from  Table  1  it  appeared  to 
give  better  overall  results.  As  in  the  previous  part  of  this  section,  a  single  gas  (N2)  was  used.  In  the 
next  section  all  work  will  concentrate  on  using  a  typical  laser  mixture.  The  calculations  reported 
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Table  3.  Five  Townsends  Comparison 

E/N 

Bin 

Max 

C020SC 

Energy 

Crompton 

(Td) 

Width(eV) 

Energy(eV) 

Balance 

Vd 

Difference 

5 

.015 

1.5 

1.136 

2.68E-05 

1.095 

3.74% 

.02 

2.0 

1.081 

3.38E-04 

1.27% 

.025 

2.5 

1.083 

4.53E-04 

1.10% 

.03 

3.0 

1.085 

5.37E-04 

0.91% 

.035 

3.5 

1.074 

1.24E-06 

1.92% 

.04 

4.0 

1.076 

3.91E-08 

1.74% 

.045 

4.5 

1.070 

4.56E-07 

2.28% 

.05 

5.0 

1.093 

8.93E-04 

0.18% 

.055 

5.5 

1.079 

8.31E^06 

1.46% 

.06 

6.0 

1.098 

1.06E-03 

0.27% 

.065 

6.5 

1.072 

1.32E-07 

2.10% 

.07 

7.0 

1.089 

3.73E-07 

0.55% 

.075 

7.5 

1.103 

1.33E-03 

0.73% 

.08 

8.0 

1.118 

2.40E-03 

2.10% 

Table  4.  Ten  Townsends  Comparison 

E/N 

Bin 

Max 

C020SC 

Energy 

Crompton 

(Td) 

Width(eV) 

Energy(eV) 

^d 

Balance 

Vd 

Difference 

10 

.015 

1.5 

2.096 

2.97E-06 

1.838 

14.04% 

.025 

2.5 

1.783 

1.35E)-04 

2.99% 

.035 

3.5 

1.785 

3.06E-06 

2.88% 

.045 

4.5 

1.786 

1.32E-06 

2.83% 

.055 

5.5 

1.805 

1.13E-06 

1.80% 

.065 

6.5 

1.802 

1.18E-07 

1.96% 

.075 

7.5 

1.840 

3.77E-04 

0.11% 

.085 

8.5 

1.811 

6.06E-07 

1.47% 

here  are  from  the  subprogram  Boltz,  using  C020SC  as  a  shell  to  run  Boltz.  This  was  accomplished 
by  setting  the  percent  of  N2  at  100,  all  other  gases  at  0,  and  using  the  nominal  settings  of  parameters 
“a”  through  “x”.  The  value  for  E/N  was  set  to  five,  ten  and  forty  Townsends  (10~^^Vcm*).  This 
was  chosen  since  good  agreement  was  achieved  over  this  range  of  E/N  in  the  previous  section  with 
other  published  results.  The  parameter  for  bin  width  was  varied  so  that  a  desired  energy  range 
was  achieved.  The  remaining  parameters  were  left  set  at  the  nominal  Values.  Table  3  lists  the  drift 
velocities  and  energy  balances  obtained. 

The  table  with  the  smallest  deviation  from  the  experimental  drift  velocity  over  all  energy 
ranges  is  Table  3,  which  is  also  the  smallest  E/N.  It  maintains  an  error  of  about  l.S%.  The  trend  is 
for  a  worsening  of  the  deviation  within  all  of  the  tables  at  the  extreme  values  of  energy  range.  The 


Table  5.  Twenty  Townsends  Comparison 


E/N 

Bin 

Max 

C020SC 

Energy 

Crompton 

(Td) 

Width(eV) 

Energy(eV) 

VJ 

Balance 

Vd 

Difference 

20 

.015 

1.5 

4.086 

1.71E-05 

3.09 

32.23% 

.025 

2.5 

3.097 

3.55E-05 

0.23% 

.035 

3.5 

3.090 

2.85E-05 

0% 

.045 

4.5 

3.107 

2.11E-06 

0.55% 

.055 

5.5 

3.135 

4.69E-07 

1.46% 

.065 

6.5 

3.148 

1.67E-07 

1.88% 

.075 

7.5 

3.197 

1.06E-07 

3.46% 

Table 

6.  Forty  Townsends  Comparison 

E/N 

Bin 

Max 

C020SC 

Energy 

Crompton 

(Td) 

Width(eV) 

Energy(eV) 

Vd 

Balance 

Vd 

Difference 

40 

.015 

1.5 

8.116 

3.52E-07 

5.18 

56.68% 

.02 

2.0 

4.640 

3.61E-06 

10.42% 

.025 

2.5 

5.285 

1.04E-05 

2.03% 

.035 

3.5 

5.295 

6.82E-09 

2.22% 

.045 

4.5 

5.342 

l.OOE-06 

3.13% 

.055 

5.5 

5.371 

2.25E-06 

3.69% 

.065 

6.5 

5.403 

1.76E-07 

4.31% 

.075 

7.5 

5.458 

3.35E-05 

5.37% 

energy  balance  is  very  good  in  all  of  the  runs.  Rockwood  (8:380)  reports  an  energy  balance  of  10~^ 
or  better.  Only  one  run  failed  to  meet  or  exceed  this  figure.  Unfortunately,  energy  balance  does 
not  appear  to  be  a  good  predictor  of  drift  velocity  computation  accuracy.  Good  estimates  and  bad 
estimates  of  drift  velocity  all  have  relatively  good  energy  balances.  Some  method  of  establishing 
the  credibility  of  a  particular  drift  velocity  estimate  is  needed.  The  motivation  for  finding  a  reliable 
indicator  is  that  it  would  help  the  user  of  C020SC  decide  if  a  particular  combination  of  energy 
bin  width  and  number  of  bins  was  a  good  choice.  This  is  important  since  drift  velocity  is  used  to 
calculate  the  total  electron  number  density,  and  this  in  turn  is  used  to  calculate  the  pumping  rates. 
In  the  cases  presented  above  the  true  answer  is  known  in  advance  and  such  an  indicator  is  not 
needed.  However,  in  those  cases  where  experimental  data  may  not  be  available  due  to  a  different 
E/N  value,  different  temperature,  or  mixture  of  gas  being  used,  such  an  indicator  would  be  very 
useful. 


A  better  indicator  appears  to  be  the  normalized  electron  number  density  distribution.  Figures 
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2  and  3  show  the  linear  and  semilog  plots  of  the  distribution  for  N2  with  E/N  equal  to  10  Td  and  a 
bin  width  of  .015  eV.  As  can  be  seen,  the  linear  plot  of  the  distribution  does  not  approach  zero  at 
the  higher  energy  bins.  It  appears  to  pass  through  a  minimum  and  start  to  rise  again.  Based  upon 
other  experiences  with  Bolts,  this  is  usually  a  good  sign  of  impending  divergence.  In  this  case  the 
estimate  of  the  electron  number  density  is  unreliable,  as  are  any  calculations  based  upon  it. 

For  the  case  where  inelastic  collisions  are  negligible  and  the  elastic  collisions  occur  with  a 
constant  collision  frequency,  a  pure  Maxwellian  would  result  (17:140).  While  it  is  expected  that 
the  inelastic  processes  present  will  drive  the  distribution  curve  away  from  a  pure  Maxwellian  form, 
the  distribution  should  still  go  to  zero  at  higher  energies.  It  is  apparent  since  the  curve  for  the 
distribution  in  Figure  2  does  not  follow  the  form  of  a  decreasing  function  at  higher  energies,  an 
energy  range  of  1.5  eV  is  inappropriate.  The  same  basic  behavior  can  be  seen  in  Figure  3. 
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In  Figures  4  and  5  are  shown  the  linear  and  semilog  plots  of  the  distribution  for  N2  wich 
E/N  equal  to  10  Td  and  a  bin  width  of  .025  eV.  The  plot  in  Figure  4  shows  a  distribution  with 
the  proper  form.  Comparison  of  the  first  two  entries  of  Table  4  show  a  marked  improvement  in  tl>e 
computational  value  of  the  drift  velocity  when  the  bin  width  of  .025  eV  is  used.  Examination  of 
Figure  5  shows  that  there  is  still  possibly  some  distribution  information  present  at  higher  energies. 

Figures  0  auid  7  show  the  linear  and  semilog  plots  for  a  bin  width  of  .035  eV  at  10  Td  E/N. 
As  can  be  seen  from  Figure  6,  it  is  impossible  to  tell  from  the  linear  plot  what  the  distribution 
function  values  are  at  energies  above  about  2  eV.  The  semilog  plot  of  Figure  7  however  reveals  tliat 
the  rate  at  which  the  distribution  decreases  has  begun  to  level  off.  The  range  of  values  for  the  plot 
in  Figure  7  covers  about  ten  decades.  Figures  8  and  9  show  a  continuation  of  this  trend  as  the 
bin  width  is  increased  to  .045  eV.  The  range  of  distribution  presented  in  Figure  9  is  about  eleven 


48 


49 


Electron  No.  Density  Nimber  of  Bins  =  108  de:  0.0358eU  h:  10.880E-08 
Linear  Plot  Tenp  :  380  Pressure  =  i  atM  E/N:  1.880H6  U  cm^2 
8.27E-81r 

8.24E-01 

8.21E-01 

0.19E-81 

8.16E-81 

0.13E-01 

8.11E-0i- 

0.88E-82 

8.S4E-82 

e.27E-02| 


9.09E*9i 


.35e£-01  .728E480  .142Et01 


.211£i01 


.281Ei01  .35eEi81 


Figure  6.  Ten  Td  Linear  Plot  for  Ae  =  .035  eV 
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decades,  and  the  drift  velocity  error  in  Table  4  also  shows  a  slight  decrease. 

The  error  tor  whe  Anal  bin  width  of  Table  4  shows  an  increase  in  error.  Examination  of  the 
linear  plot  in  Figure  10  reveals  that  the  distribution  is  squeezed  far  to  the  left  by  this  particular 
sampling  interval.  The  semilog  plot  in  Figure  11  shows  that  in  the  region  above  about  5.5  eV  the 
distribution  begins  to  exhibit  a  slight  oscillatory  behavior.  Previous  testing  with  Bolts  showed  that 
if  bin  width  was  too  large,  a  very  large  oscillatory  behavior  would  be  present  in  the  calculated  num¬ 
ber  density  distribution.  The  presence  of  this  behavior  in  Figure  1 1  is  probably  a  good  indication 
that  the  bin  width  is  too  large  and  should  not  be  used.  Table  7  presents  the  drift  velocity  error  as  a 
function  of  the  number  of  decades  from  the  distribution  maximum  to  its  minimum,  for  a  particular 
bin  width. 

While  the  data  of  Table  7  presents  some  interesting  trends  for  a  particular  E/N,  there  is 
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Table  7.  Drift  Velocity  Error  Versus  Data  Range 


E/N 

Bin 

Vd 

Decades  of 

E/N 

Bin  Vd 

Decades  of 

(Td) 

Width 

%Error 

Data 

(Td) 

Width 

%Error 

Data 

5 

.015 

3.74 

0.9 

10 

.015 

14.04 

1.0 

.02 

1.27 

3.1 

.025 

2.99 

5.6 

.025 

1.10 

10.4 

.035 

2.88 

10.2 

.03 

0.91 

17.6 

.045 

2.83 

10.9 

.035 

1.92 

20.7 

.055 

1.80 

11.2 

.04 

1.74 

21.9 

.065 

1.96 

11.3 

.045 

2.28 

22.6 

.075 

0.11 

11.4 

.05 

0.18 

26.2 

.085 

1.47 

10.8 

.055 

1.46 

22.6 

.06 

0.27 

21.9 

.065 

2.10 

21.7 

.07 

0.55 

21.6 

.075 

0.73 

21.4 

.08 

2.10 

20.9 

E/N 

Bin 

Vd 

Decades  of 

E/N 

Bin  Vd 

Decades  of 

(Td) 

Width 

%Error 

Data 

(Td) 

Width 

%Error 

Data 

20 

.015 

32.23 

0.9 

40 

.015 

56.68 

1.1 

.025 

0.23 

2.6 

.02 

10.42 

1.1 

.035 

0 

4.9 

.025 

2.03 

1.2 

.045 

0.55 

5.1 

.035 

2.22 

2.3 

.055 

1.46 

5.2 

.045 

3.13 

2.3 

.065 

1,86 

5.2 

.055 

3.69 

2.3 

.075 

3.46 

5.3 

.065 

4.31 

2.4 

.075 

5.37 

2.3 

no  one  value  with  regards  to  the  range  from  maximum  to  minimum  that  will  work  in  all  cases. 
However,  some  general  guidelines  do  appear  to  be  applicable.  First,  too  small  a  value  of  bin  width, 
resulting  in  too  small  of  an  energy  range  should  be  avoided,  or  gross  errors  in  the  calculated  number 
density  distribution  will  result.  The  energy  range  should  not  be  less  than  2.5  eV.  Likewise,  too 
large  of  a  bin  width  can  also  result  in  an  erroneous  calculated  number  density.  A  safe  upper  limit 
of  5.5  eV  appears  to  be  effective.  For  E/N,  a  range  of  from  5  to  40  Td  can  be  used.  These  limits 
also  make  sense  since  electronic  excitation  and  ionization  are  not  modeled  here.  While  all  of  these 
limits  appear  to  provide  a  good  starting  point,  the  user  of  C020SC  should  always  examine  the 
linear  and  semilog  plots  so  as  to  detect  any  signs  of  error  in  the  calculations. 
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Table  8.  Comparison  of  Predictions  for  Drift  Velocity  and  Excitation 
Elliott  Predictions  C020SC  Predictions 


E/N  (Td) 

Excitation 

Vd 

Diff 

Excitation 

Diff 

(xl0®cm/s) 

(xlO-W/s) 

(xl0®cm/s) 

(xl0~®cmVs) 

10 

3.65 

4.35 

3.58 

1.92% 

4.39 

0.92% 

20 

5.30 

5.25 

5.23 

1.32% 

5.31 

1.14% 

30 

6.70 

5.50 

6.90 

2.99% 

5.56 

1.09% 

40 

8.20 

5.60 

8.50 

3.66% 

5.71 

1.96% 

4.2  Comparison  with  Published  Results  for  a  Gas  Mixture 

In  the  previous  section,  the  transport  coefficients  calculated  by  the  Boltz  subprogram  of 
C020SC  gave  good  agreement  with  the  published  results  available  for  a  single  gas.  The  case 
of  a  typical  laser  mixture  is  considered  here.  For  this  purpose,  we  use  the  nominal  parameters 
of  C020SC.  These  settings  are  ten  percent  N2,  ten  percent  C02  and  eighty  percent  He.  The 
predictions  made  by  C020SC  are  compared  with  the  computer  predictions  published  by  Elliott 
(18:21,31),  which  were  made  with  a  program  based  in  part  on  the  work  of  Rockwood.  Comparisons 
are  made  here  between  the  computed  electron  drift  velocities  and  excitation  rates  for  the  (001) 
mode  (upper  laser  level)  of  C02.  The  results  are  listed  in  Table  8. 


We  see  in  Table  8  excellent  agreement  between  the  predictions  made  by  the  two  programs. 
We  may  now  have  confidence  that  C020SC  can  give  reasonably  good  predictions  of  the  transport 
coefficients  we  need  for  a  gas  mixture,  as  well  as  for  a  single  gas.  We  can  therefore  be  certeun  within 
the  error  bounds  listed  above  for  excitation,  that  Boltz  does  in  fact  provide  the  correct  information 
upon  which  the  pump  rate  equations  are  based. 


4.3  The  Effects  of  Boltx  Related  Parameters  on  Laser  Output 

We  have  exzunined  in  the  last  two  parts  of  this  section  the  reliability  of  Boltz  in  providing 
C020SC  the  information  needed  to  calculate  pumping  rates  baaed  upon  a  numerical  solution  of 
the  Boltzmann  equation.  This  section  examines  the  effects  on  laser  output  power  etnd  energy  of  the 
parameters  in  C020SC  that  are  directly  related  to  Bolts.  In  particular,  the  parameters  y  through 
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Table 

T(OOl)  Pulse  Length  Number 
(K)  (nsec)  of  Bins 

300  200  20 

30 

50 

100 

400  20 

30 
50 
100 

4000  200  20 

30 

50 

100 

400  20 

30 
50 
100 


.  Comparisons  of  Laser  Output 
Peak  Power  Total  Energy  density 


(MW/1) 

(J/l) 

71 

15 

71 

14 

64 

13 

64 

13 

73 

32 

72 

29 

68 

27 

69 

27 

80 

19 

81 

19 

77 

17 

76 

17 

80 

40 

81 

39 

77 

35 

81 

39 

y7  were  added  because  of  Boltz.  These  parameters  are  presented  to  the  user  in  the  opening  screen 
of  C020SC,  where  they  may  be  changed  to  whatever  values  are  needed.  We  will  now  change  the 
values  of  some  of  these  parameters  from  the  nominal  supplied  values  and  see  the  effects  as  measured 
by  the  laser  output  power  and  energy  calculations  of  C020SC. 

The  first  issue  addressed  here  is  that  of  laser  output  sensitivity  to  bin  number.  This  is  intended 
to  help  the  user  of  C020SC  to  see  how  the  number  of  bins  influences  the  predicted  performance 
of  a  particular  C02  laser  system.  Computer  runs  at  bin  numbers  20,  30,  50  and  100  were  made  at 
vibrational  temperatures  of  300K  and  4000K  for  C02  mode  (001).  This  was  done  for  pulse  lengths 
of  200  and  400  nanoseconds.  A  summary  of  the  resulting  laser  output  power  and  energy  for  each 
case  is  presented  in  Table  9. 

The  effect  of  varying  the  bin  number  (for  a  constant  energy  r2uige)  is  fairly  consistent  through¬ 
out  all  of  these  runs.  There  is  for  the  most  part  a  decrease  in  both  the  power  and  energy  as  computed 
by  C020SC  as  the  number  of  bins  increases.  The  range  of  variations  is  from  5%  to  11%  in  power 
and  from  12%  to  15%  in  energy.  The  exception  is  the  400  nsec  run  for  T(OOl)  =  4000  K.  Power  and 
energy  remain  fairly  constant  at  all  bin  numbers.  Assuming  that  as  was  previously  shown  the  100 
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bin  run  is  most  accurate,  an  approximation  of  the  error  introduced  by  using  a  smaller  number  of 
bins  is  available  from  Table  9.  If  the  answer  sought  by  the  user  is  one  which  must  predict  power  and 
energy  leveb  with  great  accuracy,  then  computer  runs  with  a  higher  number  of  bins  are  in  order. 
If  however  the  data  from  the  computer  runs  will  only  be  used  to  establish  a  trend  as  a  particular 
parameter  is  varied  over  some  range,  a  lower  bin  number  may  be  acceptrd>le. 

Examination  of  Table  9  shows  both  peak  power  and  total  energy  density  are  higher  for  the 
system  with  the  higher  effective  vibrational  temperature  for  mode  (001).  This  may  be  due  to  an 
increase  in  pump  rates,  caused  by  superelastic  collisions.  A  comparison  of  output  power  and  energy 
for  the  300  K  and  4000  K  levels  at  a  pulse  length  of  200  nanoseconds  can  be  made  by  examining 
Figures  12  through  15.  It  is  apparent  that  lasing  also  starU  sooner  for  the  4000K  system. 

The  effects  of  pulse  length  are  as  expected  for  total  energy.  The  longer  the  excitation  pulse 
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Power(HM/l)  vs  tiMe(ns)  Linear  Plot 


Figure  14.  Laser  Power  Output  with  T(OOl)  =  4000K 


Ehergy(J/l)  vs  tiMe(ns)  Linear  Plot 


Figure  15.  Laser  Energy  Output  with  T(OOl)  =  4000K 
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Table  10.  Laser  Output  Versus  E/N 


E/N 

Peak  Power 

Total  Energy  density 

(Td) 

(MW/1) 

(J/1) 

5 

14 

3.6 

10 

64 

13 

20 

75 

25 

30 

78 

28 

40 

75 

26 

lasts,  the  more  electrical  energy  will  be  converted  to  optical  energy.  That  peak  power  is  also  in¬ 
creased  when  pulse  length  is  increased  is  probably  a  result  of  gain  switching,  although  the  difference 
seen  here  is  very  small. 

The  value  of  E/N  is  seen  in  equation  (39)  and  (54a)  as  having  a  direct  influence  on  the  energy 
gained  by  the  electrons  from  the  field.  We  would  expect  from  the  cross  section  data  that  a  higher 
E/N  would  make  available  more  energy  for  transfer  from  the  electrons  to  N2.  This  can  also  be  seen 
if  the  value  of  N  is  held  constant  and  E  is  increased.  A  higher  value  for  the  applied  field  means 
a  greater  accelerating  force  is  applied  to  the  electrons.  These  more  energetic  electrons  can  then 
transfer  more  energy  to  the  other  species.  Alternatively,  if  we  hold  E  constant  and  lower  N,  the 
mean  free  path  between  collisions  will  be  increased.  This  means  that  the  electrons  will  accelerate  to 
a  higher  velocity  before  a  collision  occurs  than  in  the  case  of  the  higher  N  value.  Thus  the  electrons 
will  again  be  more  energetic  at  the  time  of  collision.  Table  10  shows  the  effects  on  laser  output  of 
varying  E/N.  The  computer  runs  of  C020SC  are  for  a  pulse  length  of  200  nanoseconds  using  100 
bins.  All  temperatures  are  at  300  K. 

The  results  presented  in  Table  10  follow  in  general  the  expected  trend.  However,  the  last  set 
of  data  for  E/N  of  40  Td  appear  to  show  a  reversal  of  the  previous  trend.  Examination  of  Figure 
16,  the  plot  of  the  normalized  electron  number  density  distribution  for  this  run,  shows  the  problem. 
The  distribution  function  is  showing  signs  of  divergence.  Switching  to  a  smaller  bin  width  would 
not  help  in  this  case,  since  the  resulting  calculation  of  the  distribution  would  not  enable  a  full 
description  of  the  entire  energy  range  covered  by  the  distribution.  This  may  also  help  explain  the 
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Figure  16.  Electron  Number  Density  for  Laser  Mixture  at  40  Td 


Table  11.  Laser  Output  Versus  Current  Density 


Current  Density 

Peak  Power 

Total  Energy  Density 

(amps/cm2) 

(MW/1) 

(J/1) 

100 

76 

17 

200 

110 

34 

observed  differences  between  the  Elliott  and  C020SC  data  of  Table  8.  Therefore  we  conclude  that 
the  40  Td  run  is  unreliable. 

The  final  Boltz  parameter  examined  here  is  current  density.  It  can  be  seen  from  section  V.A. 
that  plasma  tube  current  density  is  the  circuit  parameter  used  to  determine  the  electron  number 
density.  This  is  described  in  the  equation 

Ne  =  (current  density )/(drift  velocity  x  electron  charge) 

Thus  the  electron  number  density  is  directly  proportional  to  the  current  density.  We  would  therefore 
expect  that  as  the  current  density  is  increased,  more  electrons  will  be  present.  If  more  electrons 
are  now  available  for  collisions,  then  there  should  be  a  greater  flow  of  energy  from  the  applied  field 
to  the  excited  states  of  the  neutral  species.  This  in  turn  leads  to  the  expectation  of  greater  laser 
output  power  and  energy  when  current  is  increased,  even  when  the  other  parameters  (including 
E/N)  are  held  constant.  Table  11  shows  the  results  from  C020SC  for  two  runs  where  only  the 
tube  current  density  was  changed.  The  first  run  had  a  current  density  of  100  amps/cm^  and  the 
second  200  amps/cm^.  Both  used  100  bins  with  a  bin  width  of  .055  eV  and  E/N  of  10  Td.  Ail 
other  parameters  were  left  set  at  the  nominal  values. 

The  expected  behavior  is  seen.  More  of  the  energy  is  transferred  from  the  applied  field  to  the 
electrons  and  then  to  the  excited  states  of  the  C02  and  N2  molecules  when  plasma  tube  current 
density  is  increased.  Figures  14  and  17  show  the  peak  power  plots  for  the  100  and  200  amps/cm^ 
cases  respectively.  Aside  from  the  already  stated  difference  in  peak  power,  two  other  differences 
stand  out.  First,  the  system  with  the  higher  current  density  begins  lasing  sooner.  Since  more 
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electrons  are  available  for  impact  excitation  of  the  neutral  species,  the  inversion  needed  to  exceed 
threshold  for  lasing  can  occur  more  quickly.  The  second  difference  seen  is  in  the  backside  of  the 
main  power  peak.  The  higher  current  density  system  has  its  descent  slowed  near  the  bottom, 
resulting  in  a  somewhat  wider  pulse  than  the  lower  current  density  system.  Again,  this  is  a  result 
of  the  greater  level  of  pumping  that  takes  place  when  there  are  more  electrons  available  for  impact 
excitation. 

The  nominal  values  for  the  Boltz  portion  of  C020SC  and  suggested  operating  limitations 
are  based  in  part  on  the  results  presented  above,  and  in  part  on  the  work  of  others.  Boltz  and 
C020SC  have  no  internal  means  for  checking  on  the  reasonableness  of  the  values  a  user  might 
choose  for  the  parameters.  It  is  the  responsibility  of  the  user  to  make  sure  that  all  quantities  are 
realistic  and  not  in  conflict  with  each  other. 

The  nominal  value  for  E/N  is  10  Td.  This  value  was  tested  thoroughly  for  the  single  gas 
(N2),  and  the  mixture  C02,  N2  and  He  in  the  ratio  .1/.1/.8.  No  sign  of  divergence  or  instability 
was  observed.  The  recommended  range  for  the  gas  mixture  is  5  to  30  Td.  Above  30  Td  divergence 
appears  to  seriously  affect  the  output  for  the  gas  mixture.  It  is  not  known  if  this  is  true  for  all  other 
ratios.  The  vibrational  temperatures  of  C02  modes  (100)  and  (010)  are  set  at  the  corresponding 
ambiert  temperature  of  300  K.  The  effective  vibrational  temperature  of  mode  (001)  is  set  at  4000 
K.  Smith  (21:2042)  has  indicated  that  small  signal  gain  in  a  C02  laser  was  a  maximum  when  this 
value  vt  as  used  in  computer  simulations  he  performed.  The  number  of  bins  is  set  at  30.  The  results 
obtain  -d  indicate  that  this  is  the  best  compromise  between  speed  and  accuracy.  An  energy  range 
of  4  eV  appears  to  be  adequate  for  the  range  of  E/N  specified  above.  This  means  that  the  bin 
width  lust  be  set  to  .1333  eV.  The  nominal  value  for  current  density  comes  from  em  article  on  V-I 
characteristics  of  C02  lasers  by  Denes  (22:131). 
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V.  Conclusions 


The  subprogram  Boltz,  written  to  provide  pumping  terms  for  the  C02  laser  design  code 
C020SC,  has  been  evaluated  against  experimental  data  and  the  calculations  of  similar  programs. 
The  transport  coefficients  which  Boltz  generates  for  use  by  C020SC  appear  to  be  very  good,  when 
used  within  the  guidelines  established  in  the  previous  section.  As  noted,  the  user  of  C020SC 
must  always  be  alert  for  the  signs  of  errors  in  the  calculation  of  the  normalized  electron  number 
density  distribution.  Since  Boltz  has  not  been  tested  for  all  possible  ratios  of  C02,  N2  and  He  it  is 
impossible  to  guarantee  that  the  results  presented  here  can  be  arbitrarily  extended  to  other  cases. 
Each  individual  set  of  parameters  used  must  be  carefully  evaluated  to  prevent  an  erroneous  set  of 
data  from  causing  a  faulty  conclusion. 

The  user  of  C020SC  must  decide  the  degree  of  accuracy  with  which  the  final  laser  output 
calculations  will  be  performed.  A  more  accurate  answer  requires  more  energy  bins,  which  in  turn 
means  longer  run  times.  Two  options  allow  for  a  reduction  in  total  computer  time.  First,  initial 
runs  could  be  made  with  a  smaller  number  of  bins  in  order  to  observe  trends  in  the  laser  output 
power  and  density.  As  undesirable  combinations  are  eliminated,  a  more  accurate  computation  could 
then  be  made.  Second,  if  those  parameters  affecting  electron  kinetics  crm  be  fixed  and  run  once, 
then  parameter  y7  can  be  set  to  the  “yes”  option.  This  would  allow  subsequent  computer  runs 
to  use  the  most  recently  calculated  pumping  terms,  avoiding  the  need  for  C020SC  to  call  Boltz. 
These  two  options  can  offer  considerable  savings  in  time. 

Finally,  we  have  seen  the  effects  on  laser  output  power  and  energy  as  the  Boltz  related 
parameters  were  varied.  In  general,  results  matched  expectations  in  so  far  as  trends  were  concerned. 
No  experimental  data  and  no  calculations  from  other  programs  were  available  for  comparison  with 
the  output  from  C020SC.  In  the  case  where  an  expected  result  was  not  achieved,  the  underlying 
problem  appeared  to  stem  from  failure  of  the  calculated  distribution  to  go  to  zero  at  higher  energies. 
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A  set  of  nominal  parameters  has  been  programmed  into  C020SC.  When  used  within  the  guidelines 
stated  previously,  computational  inaccuracies  should  be  no  greater  than  those  listed  here. 
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VI.  Recommendations 


Certain  projects  could  be  undertaken  which  might  enhance  the  speed  and  accuracy  of  Boltz 
and  C020SC.  These  projects  concern  the  numerical  methods  used  by  Boltz  and  the  set  up  of 
C020SC.  Neither  Boltz  nor  C020SC  model  all  of  the  physical  processes  at  this  point  which  might 
be  of  interest. 

Among  the  methods  tested,  L-U  decomposition  offered  the  best  results  for  speed  and  accu¬ 
racy  in  solving  for  the  inverse  of  the  (I-hC)  mr.trix  However,  the  routine  used  is  general  in  its 
applicability  to  all  square  matrices.  A  Bolt;  v  •  L-U  decomposition  may  offer  a  substantial 
improvement  in  speed.  At  present,  the  L-U  decomposition  routine  may  not  be  making  best  use  of 
the  fact  that  (I-hC)  is  a  sparse,  banded  matrix.  Such  an  optimized  algorithm  may  already  exist,  or 
the  present  one  modified  to  as  to  take  advantage  of  the  special  characteristics  of  the  (I-hC)  matrix. 

The  numerical  solution  to  the  Boltzmann  equation  proposed  by  Rockwood  and  upon  which 
Boltz  is  based,  might  also  warrant  revision.  A  technique  using  an  adaptive  energy  grid  could  offer 
a  new  and  faster  approach.  This  method  offers  the  possibility  of  getting  higher  resolution  in  the 
electron  energy  distribution  in  areas  where  most  needed,  with  a  minimum  of  mesh  points  overall. 
Such  a  technique  has  been  discussed  by  Nickel  (20). 

Another  approach  which  could  be  used  to  speed  up  the  Boltz  portion  of  C020SC  would 
be  to  use  an  interpolation  table.  If  all  but  one  parameter  could  be  fixed,  a  lookup  table  could 
be  generated  over  the  dynamic  range  of  interest.  Thus,  C020SC  could  consult  the  lookup  table 
instead  of  calling  Boltz. 

Boltz  itself  does  not  take  into  account  all  of  the  physical  processes  which  a  particular  user 
might  wish  to  model.  The  effects  of  ionization  for  example  might  be  interesting  to  include.  While 
it  would  not  be  prudent  to  try  to  operate  a  laser  with  a  plasma  that  is  fully  ionized,  the  effects  of 
various  levels  of  ionization  could  be  modeled. 
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C020SC  could  also  be  improved  in  its  ability  to  model  the  effects  of  circuit  parameters.  At 
present,  only  plasma  tube  current  density  is  included.  External  excitation  by  electron  beam  would 
be  interesting  to  model.  A  technique  which  would  couple  the  external  driving  circuit  equations  to 
the  plasma  kinetic  equations  and  solve  them  in  a  seif  consistent  manner  is  one  possible  approach. 

Finally,  more  work  can  be  done  with  the  program  just  as  it  stands.  The  results  presented 
here  only  explore  the  behavior  of  C020SC  over  a  small  dynamic  range  of  the  input  parameters. 
A  more  exhaustive  set  of  test  runs  should  be  conducted  over  time  so  as  to  identify  and  document 
those  regions  where  instabilities  exist. 
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Appendix  A.  Program  Boltz 


DEFSIG  I-J 

SUB  Boltz  (kc%,  d«.  pros,  d2,  dl.  d3,  T,  VT().  El,  jd! .  Wa. 
Wb.Wc) 

OEFDBL  A-Z 

DEFSIG  D-E,  H,  J-K,  P-Q,  T.  V-W 

>*««**«  Prograa  B0LTZ31.BAS  -  coaputatioa  oi  n(u)  and 

va*^m************** 

>*«««**  Elastic,  inalaatic  and  suparalastlc  tons  for  I2,C02 

andHo  *•**•* 

*******  Incorporatos  I2DATA,  C020ATA  and  BEDATA 
files******************** 

*******  Usor  defined  Temp,  Pressnxe  ft 
E/l******************************** 


must  bebooted 
Unit. 


'lOTEM  -  arrays  are  dynaaie,  QuickBasic 
'using  qb/ak  in  order  to  exceed  64K  nenory 


PRUT 

PRIIT  "Boltznann  calculation  in  progress.” 

DIN  VB(IS) 

STATIC  etrent 

IF  event  *  1  GOTO  label4:  'test  to  see  if  this  routine 

has  already  been 

event  *  1  'run.  If  so,  do  not  re-load 

data. 

STATIC  eVO,  qij(),  eVnlO,  eVn2().  eVn3().  QlnlO,  QIn2(). 
qin3() 

STATIC  Tein>«V().  TenpOjO.  TenpyO,  jin().  I(),  ujsO, 

IsO.njsXO 

STATIC  jnl,  jn2,  jn3 

*******  Load  Cross  Section 

Data******************************************* 

OPEI  "I20ATA"  FOB  IIPUT  AS  ftl 
OPEI  "C020ATA”  FOE  IIPUT  AS  #2 
OPEI  "HEDATA”  FOR  IIPUT  AS  *3 

DIM  *V(2,  10,  SO),  qij(2,  10,  50).  *Vb1(70).  *Vb2(70). 
oVai3(70) 

DIN  qiBl(70).  qia2(70),  qia3(70)  'see 

definitions  belo« 

DIN  ToapeV(60),  Teapqj(60),  Teapy(60)  'for  use  with  linear 
interpolator 

DIN  jin(3).  K(3.  16),  uje(3.  IS),  ls<3),  BjsX(3.  IS) 

'CLS 


Road  12  inelastic  cross  section 
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data****************************** 

‘lota:  data  it  in  cn*2,  and  is  convartad  balov  to 

m*2. 

IIPUT  #1,  jin(l)  ‘jin  *  #  of  inalastic  procassas 

(axcitad  statas) 

FOR  j  =  1  TO  jin(l)  ‘nos,  stap  throngh  tha  axcitad 

statas 

IIPUT  #1.  Kd.  j).  T$,  ujsCi.  j) 

*k(l,J)  =  *  of  cross  sactions 

forthis  stata 

'Tl  -  nana  of  tha  gas  and  tha 

axcitad  stata 

‘uja  -  inalastic  thrashold 

anargy  of  stata  j 

FOR  I  =  1  T0K(1,  j)  ‘stap  throngh  tha  anargy 
lavals  of  stata  j . 

IIPUT  #1,  aV(l.  j.  I).  Qljd,  j,  I)  ‘aV  -  anargy 
inalactron  volts  for 

qijd.  j.  I)  =  Qljd,  j.  I)  /  lOOOOi  ‘anargy  I 
(orbin)  of  stata  j . 

•ext  I  *QIj  -  inalastic 

cross  saction  for 

j  ‘anargy  I  (or  bin) 

of stata  j. 

>««««««  R«ad  12  Noaantna  Transfar  Cross  Saction  Data 
a********************* 

'Iota:  data  it  in  c>*2,  and  is  convartad  balos  to 

■*2, 

IIPUT  #1,  j«l,  MTI  'j*  a  i  of  cross  tactions,  NTl  - 

nans  of  tha  gas 

FOR  j  a  1  TO  jal  ‘stap  throngh  th«  tzngj  lavals 

IIPUT  il,  aVnKj),  QlnKj)  'aV«  -  anargy  in 
alactron  volts 

qiBl(j)  a  qiBl(j)  /  lOOOOi  'qia  -  aonantna 
transfar  cross  taction 
lEZT  j 


Raad  C02  inalastic  cross  saction 


data 


‘jin  a  •  of  inalastic  procattss 


'Iota;  data  is  in  ca*2,  and  is  convartsd  balov  to 

a-2. 

IIPUT  t2,  jin(2) 

(axcitad  statas) 

FOR  j  a  1  TO  jin(2)  'now,  stap  throngh  tha  axcitad 

statas 

IIPUT  #2,  K(2,  j),  T$,  njs(2,  j) 

'k(2.J)  a  t  of  cross  sactions 

for  this  stata 

'Tl  -  nans  of  tbs  gas  and  tha 

axcitad  stata 
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'ujs  -  inelastic  threshold 


energy  of  state  j 

FOR  I  =  1  TO  K(2,  j)  'step  through  the  energy 
levels  of  state  j . 

IIPUT  #2.  eV(2,  j.  I).  qij(2.  j,  I)  'eV  -  energy 
inelectron  volts  for 

qij(2.  j.  I)  =  qij(2.  j,  I)  /  10000#  'energy  I 
(orbin)  of  state  j . 

lEXT  I  'qij  -  inelastic 

cross  section  for 

lEXT  j  'energy  I  (or  bin) 

of  state  j . 

' ******  Read  C02  MoaentUB  Transfer  Cross  Section 
Data********************* 

'lote:  data  is  in  cu*2,  and  is  converted  belos  to 

0*2. 

INPUT  #2,  jB2,  HTl  'jm  =  #  of  cross  sections,  MT$  - 

name  of  the  gas 

FOR  j  =  1  TO  ja2  'step  through  the  energy  levels 

IIPUT  #2.  eVa2(j).  qiB2(j)  'eVa  -  energy  in 

electron  volts 

qiB2(j)  *  qiB2(j)  /  10000#  'qia  •  BooentuB 
transfer  cross  section 
NEXT  j 

)««««««  R«axl  Be  NoBentua  Transfer  Cross  Section 
Data********************** 

'Bote:  data  is  in  ca*2,  and  is  converted  below  to 

a*2. 

IIPUT  #3,  ja3,  NTS  'ja  *  #  of  cross  sections,  NTS  - 

naae  of  the  gas 

FOR  j  s  1  TO  jB3  'step  through  the  energy  levels 

IIPUT  #3,  eVB3(j}.  qiB3(j)  'eVa  -  energy  in 

electron  volts 

qin3(j}  3  qin3(J)  /  10000#  'qia  -  aoaentua 
transfer  cross  section 
NEXT  j 


label4: 

FOR  I  *  1  TO  3 

input 

values 

VB(I)  ■  VT(1) 
NEXT  I 

FOR  I  *  4  TO  Jin(2) 
aabient  teap 


'set  Boltz  vibrational  teaps  to 


'set  the  reaaining  aodes  equal  to 
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VB(I)  =  T 
NEXT  I 


REDIN  abarCkc'/,).  bbar(kc'/.  1)  'for  us*  in  snergj  balance 

ft  Vd 

REDIM  A(kc'/.),  b(kcy.  +  1),  RTjs(kc*/.).  Rj8(2,  9.  kc*/.) 

REDIM  RSj8(2,  9,  kc'/.) ,  RTSjs(kc*/.)  'Suparelastics 
REDIN  c(kc'/..  kc'/.).  M0(kc'/.)  'C=  constants 

REDIN  ROl(kcX)  'NQ  =  output  energy  densities 

REDIN  IldzCkcy.),  vv(kc'/.),  yb(kc'/..  kc%)  'used  in  back-sub. 
REDIN  PlO(kcy)  'temporary  single  prec.  storage  of  ID 

for 

plotting 


'***  Defining  of  Input  Parameters  and  Constants  and  Terms  in 
S.I.  units*** 

CONST  Pi  =  3.1415926:4# 

conv  =  .00001#  'convergence  criteria  for  iteration  loop 
to  end 

'  j  =  excited  state.  s  =  species 
'  ek  s  value  of  right  hand  edge  of  energy  cell 
'Qms  -  momentum  transfer  cross  section  for  elastic 
collisions  of 

'electrons  (of  energy  ekp)  sith  molecules  of  species 
"s"  (m*2) 

'Qjs  -  inelastic  cross  section  for  this  process  (m‘2) 

'ujs  -  energy  loss  in  e-volts  for  the  j-th  inelastic  process 
in  species  s 

'T  molecular  gas  temperature  (K) 

'I  -  total  molecular  number  density  in  l/m*3 

'Nsl  -  molecular  number  density  of  species  1 

'dsl  -  number  density  of  molecules  of  species  1/1  -delta 

sub  s 

'E  =  applied  field  in  V/m 

'El  =  E/I  -  input  as  Volts  cm*2.  then  converted  to  Volts  m*2 
'de  =  width  of  cells  along  the  energy  axis  -  delta  ”*" 
(electron  volts) 

'P-  pressure  -  input  as  atmospheres,  then  converted  to 
Pascals 

q  =  (1.60220-19)  *  2  'square  of  electron  charge  -  (e*2) 
ec  3  1.6022D-19  'electron  charge 

m  =  9. 109399999999999D-31' electron  mass 
KB  =  8. 617400000000001D-06 'Boltzmann  constant  in  eV/X 
kg  =  1.3807D-23'J/X  -  S.I.  Boltzmann  constant  in  J/K 
h  *  .0000001#'tims  step 

Nsl  3  4.6618D-26  'mass  of  molecular  12  in 

Na2  =  7.3046D-26  'mass  of  molecular  C02  in  kg 

NsB  3  6.6463D-27  'mass  of  He  atom  in  kg 


P  *  pres  *  110133# 

1  =  P  /  (kg  *  T)  dcnsi'ty  in  ■*(*3) 
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EOl  s  El  /  10000* 
E  =  EOI  *  I 
dl  =  dl  /  100 
•ach  gas-  12 
d2  =  d2  /  100 
d3  =  d3  /  100 
ls(l)  =  I  «  dl 
species 

Hs(2)  =  I  e  d2 
■s(3)  =  I  «  d3 


’convert  to  V  ■‘2 
'convert  to  V/a 

'convert  to  fraction  aaount  of 

'C02 

'He 

'total  number  density  of  each 


’**•♦♦*  mjs%  =  nearest  integer  to  ujs/de  *♦♦*•♦**•*• 

'CLS 

FOR  j  =  1  TO  2  'each  species 

FOR  I  3  1  TO  jin(j)  'each  of  the  excited 

states 

mjs%(j.  I)  =  CIIT(ujs(j.  I)  /  de) 
lEXT  I 
lEXT  j 


'  *******  Calculate  Elements  of  a,  b,  and  Rjs 

Matrices**************** 

>*******  Calculate  Rja(J,K)  matrix  elements  and  the  RTjs(K) 
matrix  elements  * 

I*******  RTjs(K)  -  sum  over  all  processes 

(j) .sseeseeeeeeeseeesseeseee 

FOR  K  »  1  TP  *?  'species 

FOR  I  «  1  TO  jin(X)  'number  of  excited  states 

in  a  species 

FOR  j  =  1  TO  K(K,  I)  '*  of  data  points 
available  for  interp. 

TempeV(j)  *  eV(K,  I,  j) 

Tempqj(j)  »  qij(K,  I,  j) 
lEXT  j 

FOR  j  3  1  TO  kcX  'energy  bin 
ek  3  j  «  de 

CALL  Interp(TeapsV().  TenpQjO,  K(K.  I),  ek, 

qjo) 

RJs(K,  I,  J)  *  qjo  *  SqR(2  *  ek  *  ec  /  m>  * 

IsCK) 

RTjsCj)  3  RTjs(j)  *  qjo  *  SqR(2  *  ek  *  ee  / 

a)  *ls(K) 

CALL  Interp(TempeV() ,  TempQjO,  K(X,  I),  ek  * 
ujs(K,I).  qsjo) 

RSjs(K.  I.  j)  *  ((ek  ♦  ujs(K,  D)  /  ek)  *  qsjo  * 
SqR(2  *ek  *  ee  /  m) 

IF  K  3  1  TBEI  'selectively  apply  the  correct 

vibrational  temp 

RSja(X.  ^  j)  3  RSjs(K,  I,  j)  *  Is(K)  * 
EXP(-uj*(K.l)  /  (.  ^  T)) 

EID  IF 
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IF  K  =  2  THEI 


RSj8(K,  I,  j)  =  RSj8(K.  I,  j)  •  IsCK) 

RSj8(K,  I.  j)  =  RSj8(K.  I.  j)  ♦  EXP(-uj8(K,  I)  / 

(KB*  VB(I))) 

EID  IF 


RTSj8(j)  =  RTSj8(j)  +  RSj8(K.  I,  j) 


’*«  lota;  Rj8  and  RTj8  include  multiplication  by  la  ** 
lEXT  j 
lEXT  I 
NEXT  K 

FOR  I  =  1  TO  key. 


ak  3  I  •  da 

CALL  Intarp(aVml() ,  Qlml(),  jnl,  ak, 
CALL  Intarp(aVm2().  qim2(),  jm2.  ak. 
CALL  Intaxp(aVm3() ,  qim3(),  jm3.  ak. 
nup  3  (SqR(2  *  ac  •  ak  /  m))  *  (dl  * 
d3  *qma3) 


qmsl) 

q»«2) 

qm83) 

qmsl  *  d2  *  qm82  * 


nub  3  (<ii  a  Qaal  /  Nal)  (d2  *  qm82 
qma3  /Na3) 


nub  3  nub  «  (sqR(2  «  ac  *  ak  /  m))  « 


/  M82)  +  (d3  * 

2  *  m  *  I 


A(I)  3  (2  •  I  •  ac  /  (3  *  m))  *  {(E  /  I)  *  2)  •  (1  / 
nup)  aCda  *  -2) 

A(I)  *  A(l)  •  (ak  ♦  (da  /  4))  ♦  (nub  /  (2  •  da))  •  (KB 
•  T  /2) 

A(I)  3  A(I)  +  (nub  /  (2  *  da))  *  (-ak  ♦  (2  •  KB  •  T  / 

da)  ‘ak) 

b(I  ♦  1)  3  (2  *  I  *  ac  /  (3  ♦  m))  •  ((E  /  I)  *  2)  •  (1 


/nup)  •  (da  *  -2) 

b(I  +  1)  3  b(I  +  1)  *  (ak  -  (da  /  4))  -  (nub  /  (2  • 
da))  •(KB  *  T  /  2) 

b(I  +  1)  3  b(I  +  1)  +  (ak  ♦  (2  •  KB  •  T  /  da)  •  ak)  • 


(nub  /(2  •  da)) 
lEXT  I 


’•••••••  Load  (I-bc)-Matria  •••••••*• 

b(l)  3  0  'boundary  condition 

A(kcy.)  3  0  'boundary  condition 

c(l,  2)  3  -h  •  b(2) 

c(l,  1)  3  1  ♦  h  a  (A(l)  ♦  b(l)  ♦  RTjsd)  ♦  RTSJad)) 
c(kcy..  kcX  -  1)  3  -h  *  A(kcX  -  1) 
c(kcX.  kcX)  3  1  ♦  h  •  (A(kcX)  ♦  b(kcX)  ♦  RTj8(kcX)) 
c(kcX.  kcX)  3  c(kcX.  kcX)  *  h  *  RTSj8(kcX) 

FOR  I  3  2  TO  kcX  -  1 

c(I,  I  -  1)  3  -k  a  A(I  -  1) 
c(I,  I  +  1)  3  -h  •  b(I  +  1) 

c(I,  I)  3  1  a  h  •  (A(I)  ♦  b(I)  ♦  RTjsd)  ♦  RTSJsd)) 
lEXT  I 

FOR  K  3  1  TO  2  'spacias 

FOR  I  3  1  TO  kcX  'anargy  bin 

FOR  j  3  1  TO  jin(K)  'azeitad  state 
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IF  I  +  aj8y.(K,  j)  <=  kc%  THEM 

cim  =  -h  ♦  Rja(K,  j,  I  +  Bjsy.CK,  j)) 
c(I,  I  +  mjaJKK,  j))  *  ci«  +  c(I,  I  + 

mjsy.(K,  j)) 

EHD  IF 

IF  I  -  mjay.(K,  j)  >=  1  THEM 

ciml  =  -h  ♦  RSjsCK,  j,  I  -  ■js’/.U.  j)) 
c(I,  I  -  Bj8y,(K,  j))  =  ciMl  +  c(I,  I 

-mj8y,(K,  j)) 

EKD  IF 
lEXT  j 
EEXT  I 

MEXT  K 

'PRUT  "" 

PRUT  "  (I-h*c)  Matrix  Loadad** 


Calculata  tha  Invaraa  of 
(I-hc)*aa*******«****««******«*«« 

’***♦•♦♦*•  LU  Dacomposition  Matrix  Invaraion 
Algor itha******************* 

**********  Saa  lunar ical  Racipas  -  paga 

as******************************* 

*********  parfom  LU 

Oacoapos it ion* ************************************** 

CQIST  Tinj  *  lE-20 
0  s  1 

FOR  I  *  1  TO  kcX 
Aaaax  *  0 
FOR  j  =  1  TO  kcX 

IF  ABSCcCl,  j))  >  Aaw  THEM 
Aaaax  *  ABS(c(I,  j}) 

EID  IF 
lEXT  j 

IF  Aaaax  =  0  THEI 

PRIIT  "Singular  Matrix" 

EID 

EID  IF 

vrd)  *  1  /  Aaaax 
lEXT  I 

FOR  j  *  1  TO  kcX 

FOR  I  *  1  TO  j  -  1 
sua  »  c(I,  j) 

FOR  K  »  1  TO  I  -  i 

8ua  *  sua  -  c(I,  K)  *  c(K.  j) 
lEXT  K 

c(I,  j)  *  sua 
lEXT  I 
Aaaax  *  0 
FOR  I  »  j  TO  kcX 
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VEXT 


sua  =  c(I.  j) 

FOR  K  =  1  TO  j  -  1 

aim  =  auB  -  c(I,  K)  *  c(K,  j) 
■EXT  K 

c(I,  j)  =  auB 

duB  =  vv(I)  •  ABSCsua) 

IF  duB  >3  Aaaaz  THEM 
laaz  3  I 
Aaaaz  ~  duB 
EID  IF 
■EXT  I 

IF  j  <>  IBBX  THE! 

FOR  K  3  1  TO  kcX 

duB  3  cClBaiz,  K) 
cdaax.  K)  *  c(j ,  K) 
c( j •  K)  =  duB 
■EXT  K 
D  3  -D 

vv(Iaaz)  3  vv(j) 

EID  IF 

Ildz(j)  3  laax 
IF  c(j.  j)  3  0  THEI 
c(j.  j)  »  Tiny 
EID  IF 

IF  j  <>  kcX  THEI 

dua  3  1  /  c(j,  j) 

FOR  I  3  j  +  1  TO  kcX 

c(I,  j)  3  c(I,  j)  •  dua 
■EXT  I 
EID  IF 

j 


’ *********  Back 
Subatitution*********** 

FOR  I  3  1  TO  kcX 

FOR  j  3  1  TO  kcX 

yb(I.  j)  »  0 

■EXT  j 

yb(i.  I)  3  1 
■EXT  I 


FOR  j  3  1  TO  kcX 
II  3  0 

FOR  I  *  1  TO  kcX 
LL  3  Ildzd) 

8ua  3  yb(LL,  j) 
yb(IX,  j)  3  ybd,  j) 
IF  II  <>  0  THEI 

FOR  J1  3  II  TO  I 


'load  idantity  aatriz 


-  1 


aua  3  iUB  -  cd,  Jl)  •  yb(Jl,  j) 
■EXT  Jl 
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ELSEIF  sum  <>  0  THEI 
II  =  I 
EIO  IF 


yb(I,  j)  =  soa 
BEXT  I 

FOR  I  =  key.  TO  1  STEP  -1 
SUB  =  yb(I,  j) 

IF  I  <  key.  THEI 

FOR  J1  =  I  +  1  TO  kcX 

SUB  s  auB  -  c(I,  Jl)  *  yb(Jl, 
lEZT  Jl 
EID  IF 

yb(I,  j)  =  SUB  /  c(I,  I) 

BEXT  I 
BEXT  j 

PRIBT  "  Invars#  ol  (I-hc)  coaputad" 


j) 


’•*♦***•♦♦*  Load  initial  valuas  for  BO 


FOR  I  =  1  TO  kcX 
BO(I)  =  1 
BEXT  I 


'aaaaaaasaa  Calculata  BO  Batriz  throngk 
itarationa*aaaaa*aa*a**aaa*a***a*aa 


FOR  K  s  1  TO  300 
itarationa 


’sats  tba  saxiBUB  nuBbar  of 


FOR  I  a  1  TO  kc%  'calculata  tha  na*  alactron  nuBbar 


dansity 

SUB  a  0 

FOR  j  a  1  TO  kcX 

SUB  a  SUB  +  (yb(I,  j)  ♦  I0(j)) 
BEXT  j 

BOl(I)  a  SUB 


'PRIBT  "I  a  "  B0(I)  B0(I);  ’•  BOl(I)"; 

BOl(I) 

BEXT  I 


SUB  a  0  'calculata  tha  norBalization 

factors  for  BO  *  BOl 

SUBl  a  0 

FOR  I  a  1  TO  kcX 

SUB  3  SUB  (10(1)} 

SUBl  a  SUBl  *  (lOKD) 

BEXT  I 


cBaz  a  0  'BaziBUB  diffaraaca  valua  for  a 

particular  itaration 

FOR  I  a  1  TO  kcX  'DatarBina  coavargaaca  -  DIPT 
valuas 

nora  a  I0(I)  /  sub 

aoral  >  lOi(i)  /  subI 

diff  a  ABS(Bora  -  aoral )  /  nora 
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IF  diff  >  cmax  THEI 
cmax  =  dill 
EID  IF 

10(1)  :  101(1)  'Sat  up  10  valuas  lor  nazt  K 

iteration 
lEXT  I 

iter  =  K  ’Used  below  to  print  out  number 

ol  iterations 

IF  cmax  <  conv  THEI  'Teat  lor  convergence 

K  *  300 

EIO  IF 
lEXT  K 

FOR  I  =  1  TO  kcX 

10(1)  3  10(1)  /  aual  'normalization 

PIO(I)  3  10(1)  'storage  lor  plotting  balos.  to 
prevent  parameter 

'type  mismatch 

lEXT  I 

*******  Plot  the  normalized  electron  number  density 
distribution  *♦*♦*♦***♦ 

BEEP 

CALL  PlotKde,  kcX.  PI0(),  h,  T.  P.  EOl)  'plot  routine 
just  lor  no.  dens. 

'*••**•  Calculate  a-bar  and  b-bar  lor  drilt  velocity 
calculations  •***♦••♦• 

FOR  I  3  1  TO  kcX 
ek  3  I  «  de 

CALL  Interp(eVml(),  qimlO,  jml,  ek,  Omsl) 

CALL  Interp(eVm2() ,  QlmSO,  jn2,  ek,  Qms2) 

CALL  Interp(aVm3(),  qim3(),  ja3,  ek,  qms3) 
nup  3  (sqR(2  e  ec  «  ek  /  m))  «  (dl  *  qmsl  d2  *  qas2  '«■ 
d3  *qms3) 

nub  3  (dl  «  qmsl  /  Nsl)  *  (d2  «  qms2  /  Ns2)  *  (d3  * 
qms3  /Ns3) 

nub  3  nub  «  (SqR(2  «ec«ek/m))  *2*m*l 
abar(I)  3  (2  •  I  •  ec  /  (3  ♦  m))  •  ((E  /  I)  ‘  2)  •  (1  / 
nup)*  (de  ‘  -2) 

abar(I)  «  abar(I)  «  (ek  *  (de  /  4)) 

bbard  ♦  1)  ■  (2  *  I  *  ec  /  (3  •  m))  *  ((K  /  I)  ‘  2)  •  (1 
/nup)  •  (da  ‘  -2) 

bbard  >  1)  >  bbard  *  1)  *  (ek  -  (de  /  4}) 

IbXT  I 

'•*•••♦  Calculate  Drilt  Velocities 


CLS  0 

PRZIT  "Loop  ended  with  iteration  iter; 
PRIIT  USIIG  "  Oillerence  >  •«.•«« - cmax 
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PRIMT  •'■uab«r  of  bin*  kcX;  "  d*="; 

PRiaT  USIIG  "##.####••;  de; 

PRUT  "eV 

PRIMT  USIIG  "  I2/C02/H*  =  #.##/#.##/#.##•*;  dl;  d2;  d3 

Vd  =  0 
«bar  =  0 

FOR  I  s  1  TO  kcX  'Calculat*  nnaarical  drift 


velocity 

Vd  =  Vd  +  (abar(I)  -  bbar(I))  •  iO(I) 
ebar  =  ebar  *  (10(1)  •  (I  *  de)) 
lEXT  I 

Vd  =  Vd  *  de  /  E 

PRIMT  USIIG  "T  =####  T(IOO)  =####*•;  T;  VB(2) ; 

PRIMT  USIIG  "  T(OIO)  =####  T(OOl)  s####";  VB(1):  VB(3) 

PRIMT  USIIG  "E/M  =  . . ;  (E  /  I)  *  10000#; 

PRIMT  USIIG  "Vcm*2  Muaerical  Drift  Velocity  = 

##.### . :  Vd/.Ol; 

PRIMT  "ca/*” 

PRIMT  USIIG  "  Average  electron  energy  * 

##.#### - ebar; 

PRIMT  "eV 


' ******  Min  and  Max  Value*  of  the  Dietribution 


eeeeeeeeeeeeeeeeeeeeeeeeeeee 


Max  *  10(1) 
fory-axi* 

Min  *  10(1) 

FOR  I  *  1  TO  kcX 

IF  M0(I)  >  Max  THEM 
Max  -  M0(I) 


EID  IF 

IF  MQ(I)  <  Min  THEM 
Min  3  M0(I) 

EID  IF 
lEZT  I 


’Oeterain*  rang*  of  values 


PRIMT  USIIG  "Max  value  of  distribution  >  ««.*««« . ;  Max; 

PRIMT  USIIG  "  Min  value  of  distribution  «  #«.«#•• . ;  Min 


•«««*«•  Excitation  Rat*  Calculations  -  Fron  Rockvood  and 
Greens,  *q.(18)  ** 

•eeeeee  Mitrogen  ****** 

PRIMT  "" 

PRIMT  "  Excitation  Rates  for  Mitrogen" 

T0TI2RATE  «  0 

FOR  I  >  1  TO  jln(l)  'nunber  of  excited  states  in 

aspecies 

FOR  j  3  1  TO  K(l,  I)  ’•  of  data  points  available 
for  interp. 

T*«p*V(j)  »  *V(1.  I.  j) 

T*«pqj(j)  -  qijd,  I.  j) 

MEZT  j 
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I2IUTE  =  0 

FOR  j  =  1  TO  kcX  *  energy  bin 
ek  =  j  ♦  d« 

CALL  Int*rp(Temp*V() ,  TeapQjO,  K(l,  I),  ek,  Qjo) 
I2RATE  »  I2RATE  >  Qjo  «  SqR(2  *  ek  *  ac  /  ■)  * 

IO(j) 

HEXT  j 

PRUT  "v  =  ":  I;  ••  Rata  = 

PRUT  USIIG  ••##.#### - I2RATE  *  1000000#; 

'convaraion  to  cm*3 

PRUT  "  CB*3/8ac" 

T0TI2RATE  =  TQTI2RATE  +  12RATE  *  1000000# 

T0TI2R  s  T0TI2R  (I  *  I2RATE)  *Rata  calculation 
lor  tara  He 
lEXT  I 

PRIMT  USIHG  "Total  axcitation  rata  of  all  lavala  = 

##.### - 

T0TI2RATE; 

PRUT  "  CB-3/aac" 


C02  aaaaaa 

'Vibration  -  (010) 

FOR  j  =  1  TO  K(2,  1) 

TaBpaV(j)  =  aV(2.  1.  j) 

TaapOjCj)  s  qij(2,  1,  j) 
lEXT  j 
C02R2  »  0 
FOR  j  »  1  TO  keX 
ak  =  j  a  da 

CALL  Intarp(Ta*paV() ,  TaBpQjO,  K(2,  l),  ak,  Qjo) 
C02R2  -  C02R2  *  Qjo  *  SQR(2  *  ak  a  ac  /  ■)  *  I0(j) 
■EXT  j 

'Vibration  -  (100) 

FOR  j  =  1  TO  K(2,  2) 

TaBpaV(j)  =  aV(2,  2,  j) 

Ta*pQj(j)  =  qij(2,  2,  j) 

■EXT  j 
C02R1  -  0 
FOR  j  =  1  TO  kcX 
ak  =  j  a  da 

CALL  Intarp(TaapaV().  TaapQjO,  K(2,  2),  ak,  Qjo) 
C02R1  >  C02R1  *  Qjo  a  SQR(2  a  ak  a  .c  /  ■)  a  I0(j) 
■EXT  j 

'Vibration  -  ^wOl) 

FOR  j  »  1  TO  K(2.  3) 

Ta«paV(j)  »  aV(2,  3.  j) 

Ta«pQj(j)  .  Q*j(2.  3,  j) 

■EXT  j 
C02R3  >  0 
FOR  j  »  1  TO  kcX 
ak  a  j  a  da 
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CALL  Interp(Temp«V() ,  TempQjO,  K(2,  3),  ek,  Qjo) 

C02I13  =  C02113  +  Qjo  ♦  SQR(2  •  «k  ♦•€/■)  *  IO(j) 

MBIT  j 

PRUT  "  Ezcitatioa  rataa  for  C02" 

PRUT  USIIG  "Excitation  rata  lor  (100)  =  ##.### - C02R1 

♦1000000#; 

PRUT  "cB*3/8ac" 

PRIHT  USIIG  "Excitation  rata  lor  (010)  =  ##.### . •;  C02R2 

♦1000000#; 

PRUT  "cai“3/8ac" 

PRUT  USIIG  "Excitation  rata  lor  (001)  ♦  ##.### - ";  C02R3 

♦1000000#; 

PRUT  "cB*3/8ac" 

>«*«*««  Calculation  ol  Pumping  Tarma 
*♦♦**•**«♦«♦♦♦♦♦♦♦♦♦♦♦««♦♦♦♦♦«♦♦♦♦♦«♦♦ 

CD  =  jd!  ♦  10000!  ’Convaraion  ol  currant  danaity  to  A/m*2 
la  s  CD  /  (Vd  ♦  ac)  'Elactron  numbar  danaity  in 

l/(n-3) 

Va  ♦  la  ♦  (d2  ♦  I)  ♦  CQ2R3  'l/(8am‘3) 

Vb  a  la  ♦  (d2  ♦  I)  ♦  C02R1  'l/(8an‘3) 

He  a  la  ♦  (dl  ♦  I)  ♦  T0TI2R  'l/(aan‘3) 

>m*m***  Enargy  balaaca  chack 
♦♦♦%♦*♦♦♦♦«♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦*♦♦♦♦♦♦ 

Egain  a  o 
Ealaatic  a  0 
Einal  a  o 
FOR  I  a  1  TO  kcX 

Egain  a  Egain  ♦  (abar(I)  -  bbar(I))  •  10(1)  ♦  da 
Ealaatic  a  Ealaatic  -  (A(I)  -  abar(I)  -  b(I)  ♦  bbar(I)) 
♦  10(1)  ♦  da 

FOR  K  a  1  TO  2  'apaciaa 

FOR  J  a  1  TO  jin(K)  ’azeitad  atata 

IF  I  +  mjaXd,  j)  <a  kcX  TBEI 

Eim  a  Eja(K,  j,  I  ♦  mjaXd,  j)) 

Eim  a  Elm  ♦  io(i  ♦  mjaXd.  j))  *  njaXd. 

j)  ♦  da 

Einal  a  Einal  ♦  Eim 
EIO  IF 

IF  I  -  mjaX(K,  j)  >a  l  THEI 

Eiml  a  RSj8(K,  j,  I  -  mj8X(K.  j)) 

Eiml  a  Eiml  «  10(1  -  mj8X(X,  j))  * 

■j«X(K,  j)  ♦  da 

Einal  >  Einal  -  Eiml 
EID  IF 
lEXT  j 
lEXT  K 
lEZT  I 

PRUT  "Enargy  balanca  ia  good  to  "; 

PRIIT  USIIG  "•«.#«« - ";  ABS( (Egain  -  Ealaatic  -  Einal)  / 
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Egain) 

PRUT  "" 

IIPUT  "Press  Return  to  Continue",  11$ 

ERASE  abar,  bbar,  A.  b.  RTjs,  Rjs.  RSjs,  RTSjs,  c,  10.  101, 

Ildx,  TV,  yb,  PIO 

CLS 

EID  SUB 
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Appendix  B.  Program  C020SC  Modifications 


>mt****  C020SC  -  C02  Laser  Design  Softvare 

DECLARE  SUB  Plot2  (hi!.  h2!.  jChSt«p%,  ChTiaeStepl.  tCavl, 
jmarX, 

YquanttO.  title!) 

DECLARE  SUB  cgadump  (real) 

DECLARE  SUB  Boltz  (kcX.  del,  pres!,  d2l.  dll.  d3!.  Tl, 
VTlO. 

El!,  jd!.  Val.  Vb'.  Wc!) 

DECLARE  SUB  Interp  (z!(),  y!().  11.  xinl.  yont!) 

DECLARE  SUB  Plotl  (del,  kcX.  I0!(),  >• . .  Tl.  PI.  El!) 

DECLARE  SUB  RungeXtttta4  (tOI,  y0!().  Tl,  y!().  h!) 

DECLARE  FUlCTIOl  f!  (lY.,  T! .  y!()) 

CLEAR  ,  ,  1700  ’Alloea  increased  stack  space  for 
subprocedorea  k  functions. 

'This  C02  pulsed  laser  model  assumes  4-level  operation.  It 
incorporates 

'  choice  of  integration  time  steps  to  track  both  the 
"pulse"  and  the 

•  "tone"  portions  of  the  laser  pulse. 


>eeeee  lOMIlAL  IIPUT  PARAMETER  LIST  ***** 

Lres  s  1:  refl  s  71.6364;  areaSee  s  .0004 

pres  3  1:  pctC02  s  lO:  pctl2  =  10:  pctHe  =  80 

pctH20  3  0:  pctH2  =  0:  temp  s  300 

tauPump  3  10:  tmaz  >  400 :  hi  3  .i;  h2  3  1 

pumpEff  3  .2;  tChStep  3  30 :  ChTimeStep  «  1:  fileStep  =  5 

timeFileUnits  =  0:  isotope  3  O:  resonator  -  1:  alpha  3  3 

qSsitch  3  0:  timeqSsitch  3  15:  pulseShape  3  O:  El  3  iE-16 

’♦•**♦*  For  use  vith  Sub  Boltz 


DIM  VT(15) 

VT(2)  3  300:  VT(1)  *  300:  VT(3)  3  4OOO:  kcX  »  30:  de  3 
.1333:  jd! 

3  100 

Prevt  3  "n" 


inputroutine : 

nl  3  1;  n2  3  3  'make  array  size  variable  so  compiler 
allocates  dynamically 

IF  ChTimeStep  3  0  TBEI  n3  3  tmaz  /  hi  ELSE  n3  »  tChStep  /  hi 
♦ 

(tmaz  -  tChStep)  /  h2 
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REOIN  SHARED  Pop(nl  TO  a2,  n3) ,  Galn(n3},  Pov*r(n3), 
En«rgy(ii3) , 

REOIN  P8poiit4Pi(n3) ,  E8pont4Pi(n3)  AS  SINGLE 
REOIN  SHARED  yO(l  TO  4)  AS  SINGLE 
REOIN  SHARED  y(l  TO  4)  AS  SINGLE 
DIN  SHARED  PunpOn  AS  INTEGER 

DIN  SHARED  Ga,  Gb,  GcC02,  GcN2,  Wa,  Wb,  Wc,  Ua,  raaonatorl, 
degenRatio,  qSsitchl 

tCav  =  2  ♦  Lrea  /  cLight  /  LQG(1  /  (rafl  /  100))  'cavity 
liletima 

COLOR  15.  1 
CLS 

COLOR  16,  2 

PRINT  ••  INPUT  PARANETERS 

COLOR  IS,  1:  PRINT  " 

COLOR  14,  4:  PRINT  *'•♦•  Varaion  varaiont;  ",  Dava  Stona, 
datall;  "  ***  " 

COLOR  IS,  1 

PRINT  "a.  Raaonator  Langth  (m)  Lraa 

PRINT  "b.  Output  Railactivity  C/.)  rail;  " 

COLOR  16,  2:  PRINT  "  cavity  lilatina  = 

PRINT  USING  "«*#.»";  tCav  *  lE->'09: 

PRINT  "  nanoaac  COLOR  16,  1 

PRINT  "c.  q  Switch  on-1,  oll-O  qSwitch 

PRINT  "cl.  tiaa  to  Q  Switch  tiaaqSwitch;  "a  tCav" 

PRINT  "d.  Praaaura  (ata)  prat;  " 

COLOR  16,  2:  PRINT  "  Tha  coda 'a  tiaa  unit  =  cavity 
lilatiaa.": 

COLOR  16,  1 


PRINT  "a.  Parcant  C02 
in  V  ca*2"; 

PRINT  "  ••;  EN 

";  pctC02;  " 

y.  E/N 

PRINT  "1.  Parcant  N2 

Vibrational  Taap"; 

PRINT  "  ol  aoda  (100)";  VT(2) 

";  pctN2;  " 

yi- 

PRINT  "g.  Parcant  HE 

Vibrational  Taap"; 

PRINT  "  ol  aoda  (010)";  VT(1) 

";  pctBa;  " 

72. 

PRINT  "h.  Parcant  H2a 
Vibrational  Taap"; 

PRINT  "  ol  aoda  (001)";  VT(3) 

";  pctH20;  " 

y3. 

PRINT  "i.  Parcant  H2 

Nuabar  ol  anargy"; 

PRINT  "  bina  ";  kcX 

pctH2;  " 

PRINT  "j.  Puap  Elliciancy 

Enargy  call  width"; 

PRINT  "  in  aV  ";  da 

";  puapEll;  ■ 

•  y6. 

1 

1 

PRINT  "k.  Taaparatura  (K) 
Currant  Danaity"; 

PRINT  "  (A/ca*2)  ";  jd! 

":  taap:  " 

y6.  Tuba 
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'* :  tauPump ; 


PRI*T  "1.  Pump  Puls*  Length  (tO) 
Previous  Kinetics?"; 


"  y7. 


PRIHT  "  (y/n)  PrevI 

PRIHT  "11.  Pulse  Shape  0-rec, 1-sin  "; 
PRUT  "m.  Max  Computation  Time  (tO)"; 
PRIMT  "n.  Integration  Time  Step  #1 
PRUT  "p.  Integration  Time  Step  #2  "; 
PRIHT  "r.  Time  to  Change  Time  Step 
PRIHT  "s.  ?  change  time  step:  0,1 

PRIHT  "V.  C12  -  0  or  C13  -  1 

PRIHT  "H.  Resonator  on-1,  ofJ-0  "; 
PRIHT  "X.  Unlav  loss  y,  per  md  trip"; 
COLOR  14.  4 

IHPUT  "Select  item  to  be  changed,  <q> 


pulseShape 

tmax 

hi 

h2 

tChStep 

ChTimeStep 

isotope 

resonator 

alpha 

CO  quit,  or  <R£TURH> 


to  run:  ",  Choice! 


SELECT  CASE  Choicel 
CASE  "a",  "A" 

IHPUT  "Lres  =  ",  Lres 
CASE  "b",  "B" 

IHPUT  "refl  =  ",  refl 
CASE  "c",  "C" 

IHPUT  "qSsitch  =  ",  qSvitch 
CASE  "cl",  "Cl" 

IHPUT  "time  to  Q  Snitch  =  ",  timeqSnitch 
CASE  "d",  "D" 

IHPUT  "pres  »  ",  pres 
CASE  "e".  "E" 

IHPUT  "pctC02  =  ",  pctC02 
CASE  "f",  "F" 

IHPUT  "pctH2  =  ",  pctH2 
CASE  "g",  "G" 

IHPUT  "pctHE  =  ",  pctHe 
CASE  "h",  "H" 

IHPUT  "pctH20  =  ",  pctH20 
CASE  "i",  "I" 

IHPUT  "pctH2  »  ",  pctH2 
CASE  "j",  "J" 

IHPUT  "pnmpElf  a  ",  pumpEil 
CASE  "k",  "K" 

IHPUT  "temp  a  ",  temp 
CASE  "1",  "L" 

IHPUT  "taoPump  »  ",  taoPump 
CASE  "11".  "LI"  ' 

IHPUT  "pulse  shape  a  ",  pulseShape 
CASE  "m",  "H" 

IHPUT  "t«ax  a  ",  tmax 
CASE  "n",  "H" 

IHPUT  "hi  a  ",  hi 
CASE  "p",  "P" 

IHPUT  "h2  a  ",  h2 
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CASE  "r",  "R" 

IIPOT  "tChStep  =  ",  tChStep 
CASE  "8",  "S" 

IIPUT  "ChTia6St«p  =  ",  ChTiaeStep 
CASE  "V",  "V" 

IIPUT  "isotop*  =  ",  iaotop* 

CASE  "B",  "W" 

IIPUT  "rosonator  =  ",  resonator 
CASE  "X",  "X" 

IIPUT  "alpha  *  ",  alpha 
CASE  "y",  "Y" 

IIPUT  "E/I  in  V  cn-2  =  ",  El 
CASE  "yl",  "Yl" 

IIPUT  "T(IOO)  in  (K)  =  ",  VT(2) 

CASE  "y2",  "Y2" 

IIPUT  "T(OIO)  =  ",  VT(1) 

CASE  "y3",  "Y3" 

IIPUT  "T(OOl)  =  ",  VT(3) 

CASE  "y4",  "Y4" 

IIPUT  "luabar  oY  cells  (kcX)  »  ",  kcX 
CASE  "yS",  "YS" 

IIPUT  "call  width  de  =  ",  d* 

CASE  "y6",  "Y6" 

IIPUT  "current  density  (A/ca‘2)  =  ",  jd! 

CASE  "y7",  "Y7" 

IIPUT  "us*  previous  kinetic  calculations  (y/n)  ?  ", 

Prevt 

CASE  "q",  "Q" 

EID 
CASE  "" 

GOTO  NainBody 
CASE  ELSE 
'Continue 
EID  SELECT 

'puap  rates  and  effective  spontaneous  eaission  rat* 

d2  -  pctC02:  dl  «  pctI2:  d3  ■  pctH*  'For  us*  with  Boltz 

subprocednre 

IF  PrevI  »  "n"  THBI 

CALL  BoltzCkeX,  de,  pres,  d2,  dl,  d3,  teap,  VT(),  El, 
jd!,  Wl,  n.  W3) 

'Vb  a  Va 

>«««***  Check  to  see  if  Boltz  calculation  is  acceptable 


PRIIT  "Is  the  Boltzaann  calculation  accep4.abl*  (y/n)  ?" 
IIPUT  ”  y  -  continue  the  prograa  n  -  return  to 
aenu",  Lt 

SELECT  CASE  L$ 

CASE  "n" 

GOTO  inputroutino: 
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EID  SELECT 
ElO  IF 

PRUT  "Laser  calculation  in  progress." 

Wa  =  tCav  *  VI  /  dCrit 

Vb  =  tCav  •  W2  /  dCrit 

Vc  =  tCav  ♦  W3  /  dCrit 

Us  =  rSigna  /  areaSec  /  Pi 

Ucl  =  Uc:  Ual  =  Ua:  Vbl  =  Ub 


’Output  routines 

BEEP 

Start : 

COLOR  14,  4 

IIPUT  "Print  again?  y  or  n  ",  bt 
IF  bl  s  "y"  THEM  GOTO  Start 
titlel  =  "Poeer(MU/l)  vs  tiaeCns)" 

CALL  Plot2(hl,  h2,  jChStap,  CbTiaeStep,  tCav,  j*ai, 
PoeerO,  title!) 

titlel  =  "Energy(J/l)  vs  tiae(ns)" 

CALL  Plot2(hl,  h2,  jCbStep,  ChTiaeStep,  tCav,  jnaz, 
EnergyO,  titlel) 

IIPUT  "Make  data  files?  y  or  n  ",  cl 

IF  el  «  "y"  THEI  GOTO  OutputFiling  ELSE  GOTO  inputroutine 

RoutineVaxiableTiaeStep : 

COLOR  14,  4 
Finish: 

COLOR  14,  4 

IIPUT  "Print  again?  y  or  n  ",  bl 
IF  bl  »  "y"  THEI  GOTO  Start 
titlel  »  "Po«er(IIU/l)  vs  tiBe(ns)" 

CALL  Plot2(hl,  h2,  jChStep,  ChTiaeStep,  tCav,  jaaz, 
PoverO,  titlel) 

titlel  »  ”Energy(J/l)  tiaeCns)" 

CALL  Plot2(hl,  h2,  JChStep,  ChTiaeStep,  tCav,  jaaz, 
EnergyO,  titlel) 
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Appendix  C.  Subprograms 


DEFIIT  I-J 

DEFSIG  A-C,  F-G,  L-0,  R-S,  U,  X-Z 

CGA  Screan  Dump  Utility.  Copied  Iro*  "Microsoft 
QuickBASIC  for 

'Sciaatists"  by  Jamas  U.  Coopar,  paga  156,  Copyright  1988  by 
John  Wilay,  Inc. 

SUB  cgadump  (ras) 

DIM  buf(800)  AS  IITEGER,  I  AS  IITEGER,  j  AS  IITEGER 
DIM  col  AS  IRTEGER 
asct  s  CHRS(27} 

WIDTH  "Iptl:",  256 
SELECT  CASE  ras 
CASE  LORES: 

ctrlcharl  =  "K" 
incr  *  4 
bytas  2  400 
CASE  MEDRES: 

ctrlcharl  =  "Y" 
incr  -  8 
bytas  3  800 
CASE  HIRES: 

ctcrlchart  »  "L" 
incr  a  8 
bytas  a  800 
EID  SELECT 

LPRIIT  asci;  "3";  CHR|(24) 

FOR  col  a  0  TO  79 

DEF  SEG  a  tHBSOO 
I  a  bytas  -  1 
FOR  j  a  0  TO  99 

buf(I)  a  PEEK(j  a  80  ♦  col) 
bufd  -  1)  a  buf(I) 

IF  (ras  <>  LORES)  THEI 
bufd  -  2)  a  bufCl) 
bufCl  -  3)  a  bufd) 

EIS  IF 
I  a  I  -  incr 
lEXT  j 

DEF  SEG  a  AHBAOO 
I  a  bytas  -  1  -  incr  /  2 
FOR  j  a  0  TO  99 

bufCl)  a  PEEX(j  a  80  a  col) 
bufd  -  1)  a  bufCl) 

IF  (ras  <>  LORES)  THEI 
bufd  -  2)  a  bufd) 
bufd  -  3)  a  bufd) 

EID  IF 
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1=1-  incr 
NEXT  j 

hibyta  =  bytes  \  266 

lobyte  =  bytes  -  (hibyte  *  266) 

LPRIIT  esc$;  ctrlcharl; 

LPRIIT  CHRS(lobyta);  CHRKhibyta) ; 
FOR  I  =  0  TO  bytes  -  1 

LPRIIT  CHRKbuKD); 
lEXT  I 
LPRIIT 
MEXT  col 
LPRIIT  CHR|(12): 

END  SUB 


DEFSIG  I 

SUB  Interp  (x(),  y(),  I,  xin,  yout) 

’X  -  input  tabulated  abscissas 

(knovn) 

(knoen) 

X  *  y 
(knosn) 

(unknown) 


IF  xin  <  x(l)  THEI 
table  range, 
yout  =  0 
GOTO  label2: 
ELSEIF  xin  >  x(I)  THEI 
yout  =  0 
GOTO  label2: 

EID  IF 


’y  -  input  tabulated  ordinates 
'n  -  max  index  of  input  arrays 
'xin  -  input  abscissa  value 
'yout  -  output  ordinate  value 

'Check  to  see  if  input  x  value  is  in 
'If  it  isn't,  then  set  yout  =  0. 


klo  =  1 
khi  =  I 
labels: 

IF  khi  -  klo  >  1  THEI 

KX  >  (khi  *  klo)  /  2 
IF  x(ICX)  >  xin  THEI 
khi  >  KX 

ELSE 

klo  *  KX 
EID  IF 
GOTO  labels: 

EID  IF 

h  =  x(khl)  -  x(klo) 

IF  h  *  0  THEI 
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PRIMT  "Bad  x  input" 

EID 

EHD  IF 

yout  =  y(klo)  +  ((xin  -  x(klo))  /  h)  ♦  (y(khi)  -  y(klo)) 
lab«12: 

EIO  SUB 

DEFIIT  I-J 

SUB  Plotl  (da.  kcX.  I0().  h.  T,  P.  El) 

>*««*««  Linaar  Plot 

<«««**«  Draw  axaa  and  tick  narks 


Xmin  =  da  'Datanina  ranga  of  values 

for  x-axis 
Xmax  =  da  a  kcX 

Max  =  10(1)  'Datanina  ranga  of  values 

for  y-axis 
Min  “  0 

FOR  I  =  1  TO  kcX 

IF  10(1)  >  Max  THEI 
Max  a  10(1) 

EID  IF 
lEIT  I 

CLS  0  'Clear  the  screen 

SCREEI  2  'Sat  the  display  to  640  x  200 

UIIDOV 

VIEW 

LIIE  (64,  180) -(640,  180)  'Draw  horizontal  axis 
FOR  x  3  0  TO  10  'Draw  tick  narks  on 

horizontal  axis 

LIIE  (x  ♦  116  +  64,  l76)-(x  •  116  ♦  64,  182) 

LIIE  ((x  •  116  /  2)  ♦  64,  178)-((x  •  116  /  2)  ♦  64, 

182) 
lEZT  X 

LIIE  (64,  20)-(64,  180)  'Draw  vertical  axis 

FOR  y  3  0  TO  10  'Draw  tick  narks  on  vertical 

axis 

LIIE  (64,  180  -  y  •  16)-(68,  180  -  y  *  18) 
lEXT  y 


Label  the  Tick  Marks 


LOCATE  1.  1 

PRUT  "Electron  lo.  Density 

PRUT  "lunbar  of  bins  3";  kcX;  "  da-"; 

PRIIT  USIIG  "*«.««««";  da; 

PRIIT  "aV  h-"; 

PRIIT  USIIQ  "«##.#«# . ;  h 

PRIIT  "  Linaar  Plot  Tanp  ■  T;  "  PraSkura  ■  P  / 
110133#;  "atn"; 
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PRUT  "  E/I=": 

PRUT  USIIG  ••##.### . :  (El)  ♦  10000#; 

PRUT  '•  V  cb-2" 

YRang*  =  Max  -  Min  ’y-axis  rang# 

Ylncr  =  YRanga  /  10  '  incranaiital  diflaranca 

batsaan  y-axia  ticks 

FOR  y  =  0  TO  10  ‘Labal  vartical  axis  tick 

■arka 

LOCATE  y  *  2  ♦  3,  1 

PRUT  USIIG  "##.#####":  (Min  +  (10  -  y)  ♦  Ylncr) 
lEXT  y 
'XMin  =  0 

XRanga  >  Xaax  -  Xain  ‘x-axia  ranga 

XIncr  «  XRanga  /  S  'incraaantal  dilfaranca 

batvaan  x-axia  ticks 

LOCATE  24.  S 

FOR  X  =  0  TO  4  'Labal  horizontal  axis  tick 

marks 

PRIIT  USIIG  "««. «««««";  TAB((x)  •  14  4-  5) ;  (Xmin  *  x  * 
XIncr) ; 
lEXT  X 

PRIIT  USIIG  "  #«. #«««#“:  (Xmin  S  *  XIncr); 


Plot  tha  Data 


VIEV  (64,  20)-(639,  180)  'Oafina  via*  port  to  match  axas 
VIIDOV  (Xmin.  Hia)-(Xmax,  Max)  'Scala  pizal  eoordinatas 
FOR  X  =  1  TO  kcX  -  1 

xl  a  Xmin  *■  (XRanga  /  (kcX  -  1))  *  (x  -  1) 
yl  a  I0(x) 

x2  s  Xmin  ♦  (XRanga  /  (kcX  -  1))  *  (x) 
y2  a  I0(x  *  1) 

LIIE  (xl,  yl)-(x2,  y2) 
lEXT  X 
LOCATE  26,  1 

PRIIT  "Sand  plot  to  printar  ?"; 

LI  a  IIPUTl(l) 

SELECT  CASE  LI 

CASE  "y",  "Y" 

LOCATE  26.  1:  PRIIT  "  "; 

CALL  cgadnmp(HIRES) 

EID  SELECT 
SCREEl  0 


Log  Plot 


•aaaaaa  Dra*  axas  and  tick  marks 


Xmin  a  da  'Datarmina  ranga  ot  raluas 

for 
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XmiLX  s  d«  «  kcY. 

FOR  I  =  1  TO  kc% 

10(1)  -  .434294482#  *  LOG(IO(I))  ‘Convert  I0()  to 
common  log 
lEXT  I 


Max  =  10(1)  'Determine  range  of  values 

for  y-axis 
Min  =  10(1) 

FOR  I  =  1  TO  kcX 

IF  10(1)  >  Max  THEM 
Max  »  RO(I) 

EIO  IF 

IF  10(1)  <  Min  THE! 

Min  =  10(1) 

EID  IF 
lEXT  I 

Min  a  IIT(Min) 

Max  a  IIT(Max  ♦  1) 


CLS  0  'Clear  the  screen 

SCREEI  2  'Set  the  display  to  640  x  200 

UIIOOV 

VIEW 

LIIE  (64,  180)-(640,  180)  'Dram  horixontal  axis 
FOR  X  =  0  TO  10  'Dram  tick  marks  on 


horizontal  axis 

LIIE  (x  «  IIS  t  64. 
LIIE  ((x  *  US  /  2) 


182) 
lEXT  X 

LIIE  (64,  20)-(64,  180) 
FOR  y  a  0  TO  10 


176)-(x  •  US  64,  182) 

♦  64,  178)-((x  •  US  /  2)  ♦  64, 


'Dram  vertical  axis 
'Dram  tick  marks  on  vertical 


axis 


LIIE  (64,  180  -  y  *  16)-(68,  180  -  y  •  16) 
lEXT  y 


Label  the  Tick  Marks 


LOCATE  1,  1 

PRIIT  "Electron  lo.  Density 

PRIIT  "lumber  of  bins  a";  keX;  "  dea"; 

PRIIT  USIlO  "##.««««";  de; 

PRIIT  "eV  ha"; 

PRIIT  USIIG  "###.### - ";  h 

PRIIT  "  Semi-Log  Plot  Temp  ■  T;  "  Pressure  »  ";  P  / 
110133#;  "atm"; 

PRIIT  "  E/l«" ; 

PRIIT  USIia  "##.### . ;  (El)  •  10000#; 

PRIIT  "  V  cm*2" 

YRange  ■  Max  -  Min  'y-axis  range 

YIncr  a  YRange  /  10  'incremental  difference 
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batHeen  y-azis  ticks 

PRUT  USIIG  "»««.#«";  Max  'Labal  vartical  axis  tick 
marks 

FOR  y  =  1  TO  10 

LOCATE  y  ♦  2  +  3,  1 

PRUT  USIIG  "###.##";  (Min  +  (10  -  y)  ♦  YIncr) 
lEXT  y 
'XMin  =  0 

XRanga  »  Xmax  -  Xmin  'x-axis  rang# 

XIncr  =  XRanga  /  S  ’incramantal  dillaranca 

batsaan  x-axis  ticks 
LOCATE  24,  6 

FOR  X  =  0  TO  4  'Labal  horizontal  axis  tick 

marks 

PRUT  USIIG  ".### - TAB((x)  •  14  ♦  S) ;  (Xmin  ♦  x  • 

XIncr) ; 
lEXT  X 

PRUT  USIIG  •'  .### . ;  (Xmin  +  5  ♦  XIncr); 


Plot  tha  Data 


VIEW  (64,  20)-(639,  180)  'Oafina  via*  port  to  match  axas 
UIIDOV  (Xmin,  Min)-(Xmax,  Max)  'Seals  pixal  coordinatas 
FOR  X  *  1  TO  kcX  -  1 

xl  *  Xmin  (XRanga  /  (kcX  -  1))  *  (z  -  1) 
yl  3  |Q(z) 

x2  a  Xmin  *  (XRanga  /  (kcX  -  1))  a  (x) 
y2  =  I0(x  ♦  1) 

LIIE  (xl,  yl)-(x2,  y2) 
lEXT  X 
LOCATE  26,  1 

PRIIT  "Sand  plot  to  printar  ?"; 

Lt  s  IIPUTt(l) 

SELECT  CASE  LS 

CASE  "y",  "Y" 

LOCATE  25,  1:  PRIIT  " 

CALL  cgadnmp(HIRES) 

EID  SELECT 
SCREEI  0 
EID  SUB 

SUB  Plot2  (hi,  h2,  jChStsp,  ChTimaStap,  tCav,  jmax, 

YquantO,  titlaf) 

Linsar  Plot 


'*«««**  Dra*  axas  and  tick  marks 


GOTO  antocoords 
stplot : 

PRIIT  "Do  yon  «ish  to  spacily  th*  x  *  y  rangas  to  ba  plottad 
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?’• 

U  =  IIPUTtd) 

SELECT  CASE  L$ 

CASE  "y",  "Y" 

GOTO  mancoords 

EID  SELECT 

autocoords:  'sets  z  k  y  ranges  to  bs  plotted 

automatically 

Inin  3  0  'Determine  range  of  values 

for  z-uis 

IF  ChTimeStep  »  0  THE! 

Xmaz  =  (hi  *  tCav  *  lE't-09)  «  jmaz 
ELSE 

Xmaz  =  (jChStep  *  hi  +  (jmaz  -  jChStep)  *  h2)  ♦  tCav  • 

lE-f09 
EID  IF 

N2ut  =  Yquantd)  'Determine  range  of 

values 

for  y-azis 

Min  =  0 

FOR  I  =  1  TO  jmaz 

IF  Yquant(I)  >  Maz  THEI 
Max  s  Yquant(l) 

EID  IF 
lEXT  I 

GOTO  endcoords 
mancoords : 

IIPUT  "The  minimum  and  maximum  z  values  (ns)  are:";  Xmin, 
Xmaz 

PRUT  "Do  you  sish  to  specify  the  y  range  to  be  plotted  ?" 
LI  =  IIPUTId) 

SELECT  CASE  L$ 

CASE  "y",  "Y" 

GOTO  ycoords 

EID  SELECT 

Max  3  Vquant(l)  'Determine  range  of 

values 

for  y-axis 

Nin  3  0 

FOR  I  3  1  TO  jmax 

IF  Yquantd)  >  Max  THEI 
Max  3  Yquant(I) 

EID  IF 
lEXT  I 

GOTO  endcoords 
ycoords : 

IIPUT  "The  minimum  and  maximum  y  values  are:";  Nin,  Max 
endcoords : 

CLS  0  'Clear  the  screen 

SCREEl  2  'Set  the  display  to  640  z  200 

VIIDOV 
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'Drav  horizontal  axis 
‘Dr an  tick  narks  on 


VIEW 

LIME  (64,  180)-(640,  180) 

FOR  X  =  0  TO  10 
horizontal  aucis 

LIKE  (x  *  115  +  64.  176)-(x  ♦  116  ♦  64.  182) 

LIME  ((x  *  115  /  2)  +  64.  178)-((x  *  116  /  2)  +  64, 

182) 
lEXT  X 

LIME  (64,  20)-(64.  180)  ‘Dras  rarticad  axis 

FOR  y  =  0  TO  10  ‘Orav  tick  narks  on  vertical 

axis 

LIME  (64,  180  -  y  ♦  16)-(68,  180  -  y  ♦  16) 

MEXT  y 

Label  the  Tick  Narks 

LOCATE  1,  1 

PRUT  titlel;  "  Linear  Plot" 

YRange  =  Max  -  Min  'y-axis  range 

YIncr  =  YRange  /  10  'increnentnl  difference 

betseen  y-axis  ticks 

FOR  y  =  0  TO  10  'Label  vertical  axis  tick 

marks 

LOCATE  y  ♦  2  +  3,  1 

PRUT  USIIG  . . ;  (Min  ♦  (10  -  y)  *  YIncr) 

lEXT  y 
•XMin  »  0 

XRange  =  Xnax  -  Xnin  'x-axis  range 

XIncr  s  XRange  /  6  ‘increnentnl  difference 

between  x-axis  ticks 

LOCATE  24,  5 

FOR  X  =  0  TO  4  ‘Label  horizontal  axis  tick 

marks 

PRUT  USIIG  ".### - TAB((x)  *14+5);  (Xmin  ♦  x  • 

XIncr) ; 

NEXT  X 

PRUT  USIIG  "  .### . ;  (Xnin  +  6  *  XIncr); 


Plot  the  Data 


VIEW  (64,  20)-(639,  180)  'Define  view  port  to  natch  axes 
WIIDOW  (Xnin,  Nin)-(Xaaz,  Max)  'Scale  pixel  coordinates 
IF  ChTineStep  »  1  THEI  GOTO  varplot 
FOR  X  s  1  TO  jnax  -  1 
yl  s  Yqaant(x) 

xl  »  (hi  *  tCav  •  1E>09)  *  (x  -  1) 
x2  9  (hi  *  tCav  «  lB-t’09)  *  (x) 
y2  9  Yqnant(x  1) 

Lire  (xl,  yl)-(x2,  y2) 

lEXT  X 
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GOTO  plot end 
vairplot : 

FOR  X  =  1  TO  jmax  -  1 

IF  X  <=  jChStep  -  1  THEM 

xl  =  (hi  •  tCav  ♦  lE+09)  •  (x  -  1) 
x2  =  (hi  •  tCav  ♦  lE+09)  •  (x) 
yl  =  Yquant(x) 
y2  =  Yquant(x  +1) 

LIRE  (xl,  yl)-(x2,  y2) 

EID  IF 

IF  X  =  jChStep  THE! 

xl  =  (hi  ♦  tCav  •  lE+09)  ♦  (x  -  1) 

x2  =  (jChStep  •  hi  +  (x  -  jChStep)  •  h2)  *  tCav  * 

lE+09 

yl  ~  Yqnant(x) 
y2  ~  Yquant(x  +1) 

LIRE  (xl,  yl)-(x2,  y2) 

ERD  IF 

IF  X  >  jChStep  THER 

xl  =  (jChStep  •  hi  +  ((x  -  1)  -  jChStep)  •  h2)  • 
tCav  *  1E'*’09 

x2  =  (jChStap  •  hi  +  (x  -  jChStep)  *  h2)  *  tCav  * 

lE+09 

yl  *  Yquant(x) 
y2  s  Yquant(z  +1) 

LIRE  (xl,  yl)-(x2,  y2) 

ERD  IF 
REXT  z 
plot end: 

LOCATE  25,  1 

PRIRT  "Send  plot  to  printer 
LI  =  IRPUTt(l) 

SELECT  CASE  L$ 

CASE  "y",  "Y" 

LOCATE  25,  1:  PRIRT  " 

CALL  cgadump(BIRES) 

ERD  SELECT 

LOCATE  25,  1:  PRIRT  " 

PRIRT  "Plot  the  Linear  Plot  again  ?" 

LI  s  IRPUTKl) 

SELECT  CASE  LI 

CASE  "y",  "Y" 

GOTO  stplot 

ERD  SELECT 

Log  Plot 


•eeeeee  Drae  azee  and  tick  aarke 


FOR  I  e  1  TO  jaaz  'Convert  YqaantO  to  coaaon  log 
Yquantd)  «  .4342944821  «  LOG(Yqnant(I)} 


NEXT  I 

GOTO  autocoordsl 
stplotl : 

PRINT  "Oo  you  uish  to  spacily  tha  z  t  y  rangas  to  ba  plottad 
lor 

Log  Plot  ?" 

L$  =  INPUTSd) 

SELECT  CASE  L$ 

CASE  "y",  "Y" 

GOTO  mancoordal 

END  SELECT 

autocoordsl:  'sats  z  ft  y  rangas  to  ba  plottad 

automatically 

Xfflin  =  0  'Datarmina  ranga  of  valuas 

lor  z-azis 

IF  ChTimaStap  =  0  TEEN 

Xmaz  =  (hi  *  tCav  *  lE'fOS)  ♦  jmaz 
ELSE 

Xmaz  =  (jChStap  *  hi  +  (jmaz  -  jChStap)  ♦  h2)  *  tCav  ♦ 

lE-)-09 
END  IF 

Maz  =  Yquant(l)  'Datarmina  ranga  ol 

valuas  lor  y*azis 
Min  =  Yquant(l) 

FOR  I  =  1  TO  jmaz 

IF  Yquant(I)  >  Naz  THEN 
Max  s  Yquant(I) 

END  IF 

IF  Yquant(I)  <  Min  THEN 
Min  =  Yquantd) 

END  IF 
NEXT  I 

Min  =  INT(Min) 

Max  =  INT(Nax  *  1) 

GOTO  andcoordsl 
mancoordsl : 

INPUT  "Tha  minimum  and  maximum  x  valuas  (ns)  ara:";  Xmin, 
Xmax 

PRINT  "Do  you  wish  to  spacily  tha  y  ranga  to  ba  plottad  ?" 

L$  =  IIPUTtd) 

SELECT  CASE  Lt 

CASE  "y",  "Y" 

GOTO  ycoordsl 

END  SELECT 

Max  s  Yquant(l)  'Datarmina  ranga  ol 

valuas  lor  y-axis 
Min  3  Yquantd ) 

FOR  1  3  1  TO  jmax 

IF  Yquantd)  >  Max  THEN 
Max  3  Yquant(I) 

END  IF 
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IF  Yquantd)  <  Min  THEM 
Min  =  Yquant(I) 

EMD  IF 
NEXT  I 

Min  =  IIT(Min) 

Max  =  IIT(Max  +  1) 

GOTO  andcoordsl 


y coords 1 : 

IHPUT  "The  ninimum  and  mazimun  y  values  are:";  Min,  Max 
endcoordsl : 


CLS  0  ’Clear  the  screen 

SCREEI  2  'Set  the  display  to  640  x  200 

VIIDOW 

VIEW 

LIIE  (64,  180)-(640,  180)  'Draw  horizontal  axis 
FOR  X  =  0  TO  10  'Drav  tick  narks  on 


horizontal  axis 

LIIE  (x  ♦  116  +  64,  176)-(x  ♦  116  +  64.  182) 

LIIE  ((x  *  116  /  2)  +  64,  178)-((x  ♦  116  /  2)  ♦  64 


182) 
lEXT  X 

LIIE  (64.  20)-(64.  180) 

FOR  y  s  0  TO  10 
axis 

LIIE  (64,  180  -  y  *  16) 
lEXT  y 


'Orav  vertical  uis 
'Oraz  tick  narks  on  vertical 

(68,  180  -  y  •  16) 


•eeeeeeeeeeeeeee  Label  the  Tick  Marks 
eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee 

LOCATE  1.  1 

PRUT  title!;  "  Log  Plot" 


YRange  s  Max  -  Min 
YIncr  =  YRange  /  10 
between  y-axis  ticks 
LOCATE  3,  1 

PRUT  USIIG  "«««.»«";  Max 
FOR  y  3  1  TO  10 

narks 

LOCATE  y  *  2  3,  1 

PRIIT  USIIG 
lEXT  y 
'XNin  -  0 

XRange  »  Xnaz  -  Xnin 
XIncr  3  XRange  /  6 
between  x-axis  ticks 
LOCATE  24,  6 
FOR  X  >  0  TO  4 
narks 

PRIIT  USIIG  ".### . 

XIncr) ; 
lEXT  X 


'y-axis  range 
'increnental  difference 

'Label  vertical  axis  tick 

(Mia  ♦  (10  -  y)  •  YIncr) 

'x-axis  range 
'increnental  difference 

'Label  horizontal  axis  tick 
TAB((x)  *  14  +  6);  (Xnin  x  • 


101 


PRIHT  USIIG  "  .### . ;  (Xmin  +  S  *  XIacr) ; 


>  «**^*«***4i***«*  Plot  th6  D&ta 
««*«««««««*««««******«««««««**«««««««««««««« 

VIEW  (64,  20)-(639,  180)  'Define  viav  port  to  match  axas 
WINDOW  (Xmin,  Nin)-(Xmaz,  Max)  'Scala  pizal  coordinatas 
IF  ChTimaStap  =  1  THEN  GOTO  varplotl 
FOR  X  =  1  TO  jmauc  -  1 
yi  -  Yqnant(x) 

xl  -  (hi  ♦  tCav  ♦  lE+09)  *  (x  -  1) 

x2  =  (hi  *  tCav  ♦  lE+09)  ♦  (x) 

y2  s  Yquant(x  +  1) 

LINE  (xl,  yl)-(x2,  y2) 

NEXT  X 

GOTO  plotandl 
varplotl : 

FOR  X  =  1  TO  jmax  -  1 

IF  X  <=  jChStap  -  1  THEN 

xl  =  (hi  •  tCav  *  lE+09)  •  (x  -  1) 

x2  =  (hi  *  tCav  *  lE+09)  ♦  (x) 

yl  =  Yquant(x) 
y2  =  Yquaiit(x  +1) 

LINE  (xl,  yl)-(x2,  y2) 

END  IF 

IF  X  s  jChStap  THEN 

xl  s  (hi  *  tCav  *  lE+09)  •  (x  -  1) 

x2  *  (jChStap  •  hi  +  (x  -  jChStap)  *  h2)  •  tCav  • 

lE+09 

yl  =  Yquant(x) 
y2  =  Yquant(x  +1) 

LINE  (xl,  yl)-(x2,  y2) 

END  IF 

IF  X  >  jChStap  THEN 

xl  =  (jChStap  ♦  hi  +  ((x  -  1)  -  jChStap)  •  h2)  ♦ 
tCav  ♦  lE+09 

x2  =  (jChStap  *  hi  +  (x  -  jChStap)  *  h2)  •  tCav  • 

lE+09 

yl  =  Yqaant(x) 
y2  =  Yqaaiit(z  1) 

LINE  (xl.  yl)-(x2,  y2) 

END  IF 
NEXT  X 
plotandl: 

LOCATE  2S.  1 

PRUT  "Sand  plot  to  printar 
LI  =  IIPUTl(l) 

SELECT  CASE  LI 

CASE  "y",  "Y" 

LOCATE  26,  1:  PRUT  " 

CALL  cgadnBp(HIRES) 
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END  SELECT 

LOCATE  25.  1:  PRINT  " 

PRINT  "Plot  Log  Plot  again  ?" 
L$  =  INPUTKl) 

SELECT  CASE  L$ 

CASE  "y",  "Y" 

GOTO  stplotl 

END  SELECT 
SCREEN  0 
END  SUB 
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Appendix  D.  Data  Files 


D.l  Nitrogen 


a 

34.12  V=l,.29 
.29,0 
.3,1E-19 
.33,2E-19 
.4,3E-19 
.76,5E-19 
.9,6.5E-19 
1,8E-19 
1.1,1E-19 
1.166,1.2E-18 
1.2,1.37E-18 
1.218.1.SE-18 
1.4,6.7SE-18 
l.S,9.SE-18 
1.6,1.22E-17 
1.6S,1.39E-17 
1.7,i.eE-17 
1.8,3.3E-17 
1.9,1.52E-18 
2,1.32E-16 
2.1,4.6E-17 
2.2,1.63E-16 
2.3.1.23E>16 
2.4,4.6E-17 
2.5,8.6E-17 
2.6,1.04E-16 
2.7,2.7E-17 
2.8,4.2E-17 
2.9,4.27E-17 
3.0,4.3E-17 
3.1,6.8E-17 
3.2,3.8E-17 
3.3,2.9E-17 
3.fl,2,9E-17 
5.0,0 

18.12  Vs2,.69 
1.7,0 
1.8,9E>18 
1.9,4E-17 
2.0,1.S2E-16 
2.1,1.4E-16 
2.2.6.2E-17 
2.3,6E-17 
2.4,1.39E-16 
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2.5.1.14E-16 

2.6,3. lE-17 

2.7,4.9E-17 

2.8,5. lE-17 

2.9,1.8E-17 

3.0,2.4E-17 

3.1,1.5E-17 

3.2,1.1E-17 

3.3,7E-18 

3.4,0 

17.12  V=3,.89 

1.8,0 

1.9,1.8E-17 

2.0,7.5E-17 

2.1,1.41E-16 

2.2,1.69E-16 

2.3,8.SE-17 

2.4,2.9E-17 

2.5,7.7E-17 

2.6,1.17E-16 

2.7,6.4E-17 

2.8,2.6E-17 

2.9,4E-17 

3.0,4E-17 

3.1.1.6B-17 

3.2,1.6E-17 

3.3,1.6E-17 

3.4,0 

16.12  V*4,1.17 
1.9,0 

2.0,1.6E-17 


2.1,4.6E-17 

2.2.1.1E-ie 

2.3,1.3E-16 

2.4,7. lE-17 

2.5,2E-17 

2.6,3. lE-17 

2.7,6E-17 

2.8.4.9B-17 

2.9,1.8B-17 

3.0,1.6B-17 

3.1,1.8B-17 

3.2,1. lB-17 

3.3,7B-18 

3.4,0 

16,12  V-6.1.47 
2.0,0 
2.1,2B-17 
2.2,4.6B-17 
2.3,7.7B-17 
2.4.1.04B-16 
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2.S,1.01E-16 

2.6,6. lE-17 

2.7,2.7E-17 

2.8,3.7E-17 

2.9,6.2E-17 

3.0,4.2E-17 


3.1,2.7E-17 
3.2,3.5E-17 
3.3,3. lE-17 
3.4,0 

13,12  V=6,1.76 

2.2,0 

2.3,1.1E-17 

2.4,3.7E-17 

2.S,6E-17 

2.6,6E-17 

2.7,3.7E-17 

2.8,1.6E-17 

2.9,9E-18 

3.0,1.8E-17 


3.1,1.8E-17 

3.2,7E-18 

3.3,SE-18 

3.4,0 

12.12  V*7,2.06 
2.3,0 
2.4,7E-18 
2.S.1.8E-17 
2.6,2.9E-17 
2.7,4.4E-17 
2.8,3.3E-17 
2.9,1.8E-17 
3.0,6E-18 
3.1,7E-18 
3.2,1.6E-17 
3,3,7E-18 
3.4,0 

8.12  Vs8,2.36 
2.6,0 
2.6,7E-18 
2.7,1. lB-17 
2.8,1.8E-17 
2.9,2.48-17 
3.0,1.68-17 
3.1,78-18 
3.2,0 

60.12 
0,18-16 
18-4,18-16 
48-4.1.28-16 
98-4,1.338-16 
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1.6E-3,1.43E-16 

2.5E-3,1.56E-16 

3.6E-3.1.69E-16 

4.9E-3.1.8E-16 

6.4E-3,1.94E-16 

8.1E-3,2.05E-16 

1.03E-2.2.2E-16 

1.44E-2,2.49E-16 

2.21E-2,2.94E-16 

3.24E-2,3.SE-16 

4.0E-2,3.8«E-ie 

4.84E-2.4.24E-16 

6.51E-2.4.9E>16 

7.86E-2,5.33E-16 

.103,6.04E-16 

.11S6.6.31E-16 

. 1502,7. 12E-18 

.226.8.22E-16 

.332,9.34E-16 

.446.9.9SE-16 

1.0.9.98B-16 

1.1.1.014E-1S 

1.2.1.061B-1S 

1.3,1.18E>16 

1.4,1. 146E-16 

l.S,1.196E-16 

1.6,1.29E-1S 

1.7,1.343E-15 

1.8.1.696E-1S 

1.9.1.983E-16 

2.0,2.401E-15 

2.2,2.876E-15 

2.6,2.9888-16 

2.8,2.801E-1S 

3.0,2. 163E-16 

3.3,1.719E-16 

3.6,1.466E-16 

4.0,1.262E-16 

4.6,1. 162E-15 

6.0.1.03E-1S 

10,8.51E-16 

15,1.1E-15 

20,1.2E-16 

26,1.17E-16 

36,1.06E-16 

40,1.01E-16 
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D.2  C02 


9 

36.C02  vs 1. .0832 

0.0 

.0827,0 
.0844,8.5E-17 
.0862.1.16E-16 
.0932.1.8SE-16 
.1035.2.3E-16 
.1208,2.6E-16 
.1382.2.68E-16 
.1726.2.62E-16 
.207,2.48E-16 
.275,2. 18E- 16 
.346,1.93E-16 
.S.1.45E-16 
.7.1.1E-16 
.9.8E-17 
1.1.6.2E-17 
1.4,4.6E-17 
1.6,4.2E-17 
1.8,4.4E-17 
2.3,7E-17 
2.6,9.3E-17 
3.0,1.34B>16 
3.2.1.S8E-16 
3.4,1.7SE-16 
3.6.1.8E-16 
3.8.1.79E-16 
4.0,1,7E-16 
4.2,1.62E-16 
4.6.1.06E-16 
S.1,5.7E-17 
5.5,5. lE-17 
6,5E-17 
7.4.8E-17 
8,4.5E-17 
10,2E-17 
20,0 

27.C02  v»2,.167 

0.0 

.167,0 

.2,6E-17 

.22,7.6E-17 

.26,8B-17 

.3,7.8B-17 

.5,6.4B-17 

.7.5.3E-17 

1.4.4E-17 
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4.5,1.6E-17 

S.06.4.4E-18 

6.0,0 

5,C02  v=6,.422 
2.5,0 

3.0,1.0SE-17 

3.S6,2.25E-17 

4.1,1E-17 

4.5,0 

5,C02  v=7,.505 
2.5,0 

3.0,1.56E-17 

3.56,3.3E-17 

4.1,1.56E-17 

4.5,0 

5,C02  v==8,2.5 
2.5,0 

3.0,1.8E-17 

3.8,2.5E-17 

4.1,i.8E-17 

4.5,0 

10,C02  v=e,3.85 
3.65,0 
4.3,1.4E-19 
4.5,1.4E-19 
5.1,0 
6.6,0 


7.2,7E-20 

8.2,4.6E-19 

8.4,4.2E-19 

8.9,1E-19 

9.7,0 

44,C02 

0.0,6E-14 

.001,5.4E-14 

.002,3.8E-14 

.004,2.7E-14 

.007,2E-14 

.01,1.7E-14 

.02.1.2E-14 

.04,8.6E-15 

.07,e,4E-18 

.1,5.2E-16 

.15,4E-16 

.2,3.16E-16 

.26,2.5E-15 

.3,2E-1S 

.36.1.66E-16 

.42,1.3E-16 

.5,1.08E-16 

.6,8.8E-16 


no 


.7,1.1E-16 

.85.6.3E-16 

1.0,S.6E-16 

1.25,6. lE-16 

1.5,SE-16 

1.8,6E-16 

2.2,S.4B-16 

2.5,6.5E-16 

2.8.7.6E-16 

3.2,1.03E-15 

3.6,1.42E-15 

4.0,1.52E-1S 

4.5,1.48E-15 

4.9,1.32E-16 

5.2.1.2E-1S 

S.6.1.08E-1S 

6.4,9. 8E-16 

8.0,1.08E-16 

10,1.21E-16 

14.1.41E-1S 

18,1.66E-1S 

28,1.62E-1S 

40.1.46E-1S 

S2,1.27E-16 

7S,9.6E>16 

100,8E-16 


D.3  Helium 


38.  H« 

0.0,5. 18E-16 

.008,5. 18E- 16 

.009,5.198-16 

.01.S.21E-16 

.013.6.2eE-ie 

.017,5.31E-16 

.02,6.36B-16 

.025.5.41E-16 

.03.6.46B-16 

.04.6.64B-16 

.06,6.62B-16 

.06.6.6BB-16 

.07,6.74B-16 

.08,5.79B-17 

.09.6.83E-ie 

.1,6.86E-16 

.12,6.94B-16 

.16.6.04B-16 
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Appendix  E.  Methods  of  Numerical  Solution 


E.l  Gauss- Jordan 


The  Gauss-Jordan  method  calculates  the  inverse  of  a  matrix  through  a  series  of  steps  called 
normalization  and  reduction  (11:193).  The  first  step  in  this  process  is  to  form  an  augmented  matrix 
using  the  A  matrix  which  is  to  be  inverted,  and  matrix  /,  the  identity  matrix. 


ail 

012 

013 

1 

0 

0 

021 

022 

023 

0 

1 

0 

(51) 

asi 

032 

O33 

0 

0 

1  ; 

The  second  step  is  to  normalize  all  the  elements  of  the  first  row  by  its  diagonal  element. 
The  normalized  first  row  is  now  multiplied  by  an  appropriately  chosen  constant  such  that  when 
subtracted  from  the  second  row  the  result  will  reduce  at  least  one  of  the  second  row’s  elements 
to  zero.  This  process  of  row  reduction  using  the  normalized  first  row  is  performed  on  each  row 
(excluding  the  normalized  row).  When  reduction  using  the  first  row  is  complete,  the  second  row  is 
normalized  by  its  diagonal  element,  and  reduction  of  all  the  other  rows  performed.  This  process  is 
repeated  until  the  last  row  has  been  normalized.  Once  complete,  the  original  matrix  A  will  have 
been  reduced  to  the  identity  matrix,  and  the  /  matrix  will  have  been  replaced  by  the  inverse  of  A. 
This  method  can  be  used  to  solve  for  the  inverse  of  (/  -  hC)  in  equation  (42).  An  initial  guess 
distribution  is  then  used  to  iteratively  solve  for  the  electron  number  density  distribution.  The  initial 
guess  distribution  used  was  a  straight  line  distribution  (constant  value  equal  to  one).  Fortunately, 
this  inverse  need  only  be  calculated  once  for  a  particular  solution.  The  iterative  process  of  equation 
(42)  is  repeated  until  the  desired  convergence  criteria  (10~’)  is  met.  The  algorithm  used  for  the 
Gauss-Jordan  solution  was  taken  from  reference  11. 
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E.2  L-U  Decomposition 


L-U  decomposition  is  another  method  which  attempts  to  solve  equation  (42)  by  calculation 
of  the  inverse  of  (7  —  hC).  This  technique  is  based  upon  a  decomposition  of  the  original  matrix  A 
into  a  lower  triangular  part  and  an  upper  triangular  part.  These  are  denoted  L  and  U ,  such  that 

LU  =  A  (52) 

and  we  wish  to  solve  the  linear  equation 

Ax  =  b  (53) 

such  that 

i  f  =  (f  t/)  •£=  L- (^  x)  =?  (54) 

We  accomplish  this  by  first  solving  for  the  vector  y  where 

Ly  =  b  (55) 

and  then  solving  for  x  by 

Ox  =  y  (56) 

Once  we  have  L-U  decomposed  a  particular  matrix  A,  we  can  solve  for  its  inverse  one  column 
at  a  time,  by  successively  letting  the  6  vector  be  a  column  of  the  identity  matrix.  The  resulting 
columns  when  put  together  will  form  the  inverse  of  the  A  matrix.  The  program  for  this  procedure 
comes  from  Press  (10:31).  This  inverted  (I  —  hC)  matrix  was  then  used  iteratively  in  equation  (42) 
with  the  initial  guess  being  a  straight  line  (constant  distribution). 
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E.3  Gauss-Seidel 


Referring  to  equation  (40),  we  see  that  for  the  steady  state  case 


^  Ctin,  =  0  (57) 

t 

Assuming  once  again  that  the  elements  of  the  C  matrix  are  constants,  we  are  left  with  solving  a 
set  of  linear  equations.  Iterative  algorithms  can  offer  a  fast,  accurate  solution  in  cases  where  the  C 
matrix  is  large  (13:382).  One  such  technique  is  the  Gauss-Seidel  algorithm.  For  solving  equation 
(49),  this  algorithm  takes  the  form  (9:373) 


^  fr  +  E  ■’')  (58) 

*'  Vj=1  /=*+!  / 

foT  k  =  K  to  k  =  1.  This  is  an  improvement  over  a  technique  such  as  the  Jacobi  method  because  as 
soon  as  they  are  available,  updated  values  are  used  to  calculate  in  the  second  term  on  the  right  to 
calculate  estimates  of  the  same  iterative  order.  The  rate  of  convergence  may  be  improved  through 
the  use  of  Aitken’s  acceleration  technique  (9:373).  The  formula  for  this  algorithm  is 


„«+i)  =  „(«-2) 


(59) 


A  potential  draw  back  with  the  Gauss-Seidel  method  lies  with  the  initial  guess  that  must  be 
supplied.  If  the  initial  guess  is  far  off,  convergence  may  be  very  time  consuming  (even  with  the 
use  of  an  accelerator).  In  order  to  limit  this  problem,  the  input  initial  guess  used  came  from  20, 
30  and  50  bin  calculations  using  the  L  U  decomposition  method.  An  initial  guess  constant  vrdue 
distribution  was  also  tried. 
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E.4  Successive  Over  Relaxation  (SOR) 


Let  us  denote  x  as  an  approximation  to  a  linear  system  Ax  =  6.  VVe  will  define  a  residual 
vector  with  respect  to  this  system  as  f  =  b  —  Ax.  Using  the  notation  of  the  residual  vector,  the 
Gauss-Seidel  method  may  be  rewritten  as  (13:391) 


(‘fc) 

an 


This  in  turn  may  be  modified  such  that 


Certain  choices  of  w  may  lead  to  faster  convergence  for  equation  (53).  Methods  of  solution 
based  upon  equation  (53)  are  called  relaxation  methods.  Under-relaxation  is  the  region  where 
0  <  w  <  1.  Over-relaxation  is  the  region  where  w  >  1.  Both  of  these  methods  are  abbreviated  as 
Successive  Over- Relaxation  (SOR).  Using  this  approach,  we  rewrite  equation  (50)  as  (13:391) 


>=i  ;=i+i 


for  i  =  1  to  K.  The  advantage  of  this  method  over  Gauss-Seidel  is  that  it  may  provide  accelerated 
convergence.  According  to  the  theorem  of  Kahan  (13:392),  we  must  have  0  <  u;  <  2  if  Oi,  ^  0  for 
each  I  =  1,2,  ...K.  Otherwise  the  solution  will  not  be  convergent.  As  in  the  case  of  the  Gauss- 
Seidel  method,  L-U  decomposition  and  a  straight  line  distribution  were  used  to  seed  the  initial 
guess  vectors. 
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E.5  Speed  and  Accuracy  Comparison 

To  judge  how  well  a  particular  method  worked,  a  normalized  n(e)Af  =  1)  electron  number 
density  distribution  was  calculated  using  the  above  methods.  Double  precision  was  used  in  all 
calculations.  The  computer  used  was  an  Apple  IIGS  with  an  Applied  Engineering  PC  Transporter, 
math  coprocessor  and  640K  of  ram  installed.  The  cross  section  data  used  was  for  nitrogen  (14:62- 
67).  From  this  distribution  function,  the  drift  velocity,  relative  convergence  of  the  distribution 
and  the  energy  balance  were  calculated.  To  determine  the  energy  balance,  the  rate  at  which  the 
electrons  gained  energy  from  the  field  and  the  rates  at  which  the  electrons  lost  energy  through 
elastic  and  inelastic  collisions  were  all  calculated.  From  equation  (39),  (56)  and  (57)  we  get 

Ef  =  -  hk)nic^( 

k 

Eel  =  —  -  ait)  -  {hk  -  it)]nk  Ac 

k 

Eint!  —  ~  ~  ^jtE!j,nit  —  m^'f)]Ac 

k 

Energy  Balance  =  ^ ^ (63) 

A  convergence  value  for  each  element  of  the  distribution  was  calculated  according  to 
|nj*  -  The  convergence  criteria  used  was  10"®. 

The  number  of  iterations  and  run  time  were  also  recorded.  The  calculated  drift  velocity  was 
compared  with  that  for  an  experimentally  reported  value  (15:627-628)  with  E/N  =  5x10"^^  V  cm^ 
at  a  temperature  of  293  K.  The  results  are  reported  in  Table  12. 

In  Table  12,  the  best  results  obtained  by  each  method  are  reported.  The  two  methods  used 
to  seed  the  initial  guess  for  the  last  three  methods  were  L-U  decomposition  and  a  straight  line 
distribution.  Runs  were  made  using  20,  30  and  50  bins  to  generate  the  initial  guess  distribution 
by  L-U  decomposition.  The  best  results  were  obtained  using  fifty  bins,  although  this  increased  run 
time.  L-U  decomposition  providing  the  initial  guess  always  gave  better  speed  and  equal  or  better 
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Method 

vd*** 

Table  12.  Comparison  of  Algorithms 
Converg  Balance  Iter  Run  Time 

Experimental 
L-U  decomp 

1.095 

1.083 

1.28x10“^ 

4.50xl0~‘‘ 

4 

9.5 

Gauss-Jordan 

1.083 

2.76x10“® 

4.50x10“* 

4 

18 

Gauss-Seidel* 

1.079 

9.98x10“® 

4.51x10“® 

1188 

90 

Gauss- Seidel** 

— 

D 

— 

— 

— 

SOR 

1.081 

9.86x10“® 

1.48x10“* 

112 

10.5 

*  -  no  acceleration  technique  is  applied 
**  -  Aitken’s  acceleration  technique  is  applied 
***  -  velocity  is  xlO®  cm/s 
D  -  did  not  converge 

All  run  times  are  in  minutes,  convergence  criteria  was  10“  ® 


accuracy  than  the  straight  line  distribution. 

Gauss-Seidel  with  Aitken’s  acceleration  did  not  converge  to  a  stable  solution.  An  oscillatory 
component  was  present  in  the  convergence  calculations,  and  this  prevented  the  calculation  from 
establishing  a  sufficiently  strong  overall  downward  trend.  This  same  result  occurred  regardless  of 
the  initial  guess  used. 

SOR  also  exhibited  an  oscillatory  behavior  in  the  convergence  calculation.  However,  this 
behavior  was  weak  enough  so  that  the  overall  trend  in  convergence  was  downward.  Eventually,  the 
convergence  criteria  was  met.  The  best  SOR  results  were  obtained  with  a  50  bin  L-U  decomposition 
initial  guess  and  cv  =  1.9. 

L-U  decomposition  and  Gauss-Jordan  provided  results  of  comparable  accuracy.  Run  time  for 
Gauss-Jordan  is  however  nearly  twice  that  of  L-U  decomposition.  This  factor  of  two  is  due  to  the 
greater  number  of  computations  required  of  the  Gauss-Jordan  method.  About  eighty  five  percent  of 
run  time  was  spent  on  matrix  inversion  for  both  L-U  decomposition  and  Gauss-Jordan.  Experience 
has  shown  that  an  energy  balance  of  10“^  or  smaller  (meaning  a  better  energy  balance)  should  be 
obtainable  for  a  good  choice  of  energy  range,  bin  width,  number  of  bins  and  time  step.  An  energy 
balance  of  greater  than  10“^  has  usually  meant  the  presence  of  an  instability  which  might  be  seen 
by  examination  of  the  linear  and  log  plots  of  the  number  density.  Of  all  the  methods  tried,  L-U 
decomposition  provided  the  best  accuracy  and  run  time. 
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